clear all
global a b 

a=0.01;

Tfin=200*pi;
e=0.01;
ei2=1/sqrt(5)*[1;2];
ei1=[1;0];
eid=[0;1];
x0=4;y0=0;

dt=0.01;
Nt=ceil(abs(Tfin/dt));
xx=zeros(Nt,1);yy=xx;tt=xx;
xx(1)=x0;yy(1)=y0;tt(1)=0;
xt=x0;yt=y0;

for it=2:Nt;
    tt(it)=tt(it-1)+dt;
[k1x,k1y]=synartisi2d(xt,yt,tt);
[k2x,k2y]=synartisi2d(xt+k1x*dt/2,yt+k1y*dt/2,tt+dt/2);
[k3x,k3y]=synartisi2d(xt+k2x*dt/2,yt+k2y*dt/2,tt+dt/2);
[k4x,k4y]=synartisi2d(xt+k3x*dt,yt+k3y*dt,tt+dt);
xt=xt+dt*(k1x+2*k2x+2*k3x+k4x)/6;
yt=yt+dt*(k1y+2*k2y+2*k3y+k4y)/6;
xx(it)=xt;yy(it)=yt;
end;

figure(10);
plot(xx,yy,'r');
axis('square')
hold on