clear all;


figure(1);
hold on;

x0=1;
y0=1;



tf=25;
Nt=801;

T=linspace(0,tf,Nt);
d= T(2)-T(1);




X=zeros(2,length(T));


X(:,1)=[x0;y0];


for k=2:length(T);
    
    t=T(k-1);
    
    x=X(1,k-1);
    y=X(2,k-1);
    
    [ky1,kv1] = twoD( x         , y         , t     );
    [ky2,kv2] = twoD( x+d*ky1/2 , y+d*kv1/2 , t+d/2 );
    [ky3,kv3] = twoD( x+d*ky2/2 , y+d*kv2/2 , t+d/2 );
    [ky4,kv4] = twoD( x+d*ky3   , y+d*kv3   , t+d   );
    
    
    xnew = x + d/6*(ky1+2*ky2+2*ky3+ky4);
    ynew = y + d/6*(kv1+2*kv2+2*kv3+kv4);

    X(:,k)=[xnew;ynew];
     
end


hold on;
plot(X(1,1:k),X(2,1:k),'b')
% hold on;
% plot(X(1,1),X(2,1),'.r','MarkerSize',8)
hold off;
axis square;
axis([-4 4 -4 4]*2)
xlabel('x','FontSize',18);
ylabel('y','FontSize',18);
title('phase space','FontSize',18);  
drawnow;