clear all;

N = 1;
h = 2*pi;




x = [1];
y = [0];
xt = 0;



hold on;

th = 0:pi/50:2*pi;
xunit = cos(th);
yunit = sin(th);
plot(xunit, yunit);

for i=0:N*101
    xt=x(end);
    yt=y(end);
    x = [x,(xt-h*yt-1/4*h*h*xt)/(1+1/4*h*h)];
    y = [y,1/2*h*(x(end)+xt)+yt];
end

axis([-1.1 1.1 -1.1 1.1]);
axis square;
plot(x,y,'r');

hold off;
