clear all
 
global epsilon

epsilon=2;

a=.1;

dt=0.001;

Nt=100000;

psi=[0;2];

tt(1)=0;x(1)=psi(1);y(1)=psi(2);



for it=2:Nt
    tt(it)=(it-1)*dt;
    t0=tt(it);
    k1=dt*xdot2G(psi,t0);
    k2=dt*xdot2G(psi+k1/2,t0+dt/2);
    k3=dt*xdot2G(psi+k2/2,t0+dt/2);
    k4=dt*xdot2G(psi+k3,t0+dt);
    psi=psi+(1/6)*(k1+2*k2+2*k3+k4);
    %psi=psi+k1;

    x(it)=psi(1);y(it)=psi(2);

end
%%
clear tt x y

psi0=psi;

tt(1)=0;x(1)=psi(1);y(1)=psi(2);
PHI=eye(2);

for it=2:Nt
    tt(it)=(it-1)*dt;
    t0=tt(it);
    k1=dt*xdot2G(psi,t0);
    k2=dt*xdot2G(psi+k1/2,t0+dt/2);
    k3=dt*xdot2G(psi+k2/2,t0+dt/2);
    k4=dt*xdot2G(psi+k3,t0+dt);
    psi=psi+(1/6)*(k1+2*k2+2*k3+k4);
    %psi=psi+k1;

    x(it)=psi(1);y(it)=psi(2);
    et=norm(psi-psi0);
    J=[0,1;-1,epsilon*(1-(y(it))^2)];
    AE=expm(dt*J);
    PHI=AE*PHI;


    if et<1.e-3
        break
    end
end

T=tt(it);

T/(2*pi)

mu=abs(eig(PHI));

lambda=log(mu)/T
lambda/epsilon



ta=linspace(0,2*pi,401);
xa=2*sin(tt);ya=2*cos(tt);
figure(1)
plot(tt,x,'r',tt,xa,'--b','linewidth',2)
xlabel('$t$','Interpreter','latex')
ylabel('$x(t)$','Interpreter','latex')
xlim([tt(end)-4*pi,tt(end)])
set(gca,'Fontsize',24)

figure(2)
plot(x,y,'r',xa,ya,'--b','linewidth',2)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
set(gca,'Fontsize',24)




    
