%program for integration 2D dynamical systems

global e a

e=0.02; a=1;


xt0=[0;3];

dt=0.1;

Nt=100*ceil(2*pi/dt);

xt=zeros(2,Nt);
tt=zeros(Nt,1);

xt(:,1)=xt0;

for it=2:Nt;
tt(it)=(it-1)*dt;

xtp=xt(:,it-1);ttp=tt(it-1);

k1=dt*fdot2D(xtp,ttp);
k2=dt*fdot2D(xtp+k1/2,ttp+dt/2);
k3=dt*fdot2D(xtp+k2/2,ttp+dt/2);
k4=dt*fdot2D(xtp+k3,ttp+dt);

xt(:,it)=xtp+(1/6)*(k1+2*k2+2*k3+k4);
if rem(it,100)==0
figure(33)
plot(xt(1,1:it),xt(2,1:it),'r','linewidth',2)
xlabel('$x$','interpreter','latex')
ylabel('$y$','interpreter','latex')
set(gca,'Fontsize',18)
axis('equal')
drawnow
end

end;

figure(35);hold on
plot(xt(1,1:it),xt(2,1:it),'b','linewidth',2)
xlabel('$x$','interpreter','latex')
ylabel('$y$','interpreter','latex')
set(gca,'Fontsize',18)
axis('equal')
drawnow



% figure(34)
% plot(tt,xt(1,:),'r',tt,exp(e/2*tt).*cos((sqrt(1-e^2/4))*tt),'.b','linewidth',2)
% xlabel('$x$','interpreter','latex')
% ylabel('$y$','interpreter','latex')
% set(gca,'Fontsize',18)
% drawnow

