%program for integration 2D dynamical systems
clear all
global mu nu  xe1 xe2

%c <1
% nu^2 > 4 mu 

c=1/2;

mu=1/2;
nu=1/4;
%nu=sqrt(3);

[nu^2,4*mu]
pause
K=(1-c^2)/2;
xe1=c-1
xe2=c+1

s=1;
A=[0,1;s/mu,nu/mu];
[u1,s1]=eig(A)
una=u1(:,2) %unstable
sta=u1(:,1) %stable
% s=-1;
% A=[0,1;s/mu,nu/mu];
% [u2,s2]=eig(A)


dt=-0.001;


Nt=19*10500;


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

xt0=[xe1;0]-0.000001*u1(:,1);


xt(:,1)=xt0;

for it=2:Nt

tt(it)=(it-1)*dt;
tp=tt(it);
xtp=xt(:,it-1);

k1=dt*fdot2DKdVB(xtp,tp);
k2=dt*fdot2DKdVB(xtp+k1/2,tp+dt/2);
k3=dt*fdot2DKdVB(xtp+k2/2,tp+dt/2);
k4=dt*fdot2DKdVB(xtp+k3,tp+dt);

xt(:,it)=xtp+(1/6)*(k1+2*k2+2*k3+k4);

end








figure(306);hold on
plot(xt(1,1:it),xt(2,1:it),'r','linewidth',2)
xlabel('$u$','interpreter','latex')
ylabel('$\dot u$','interpreter','latex')
set(gca,'Fontsize',24)
%axis('equal')
box on
drawnow


figure(1);hold on
plot(tt,xt(1,:),'r','linewidth',2)
xlabel('$\xi$','interpreter','latex')
ylabel('$u(\xi)$','interpreter','latex')
set(gca,'Fontsize',24)
%axis('equal')
box on
xlim([-50 0])
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
% figure(34)
% plot(tt,xt(1,:),'r',tt,sqrt(12)*cos(tt),'b','linewidth',2)
% xlabel('$t$','interpreter','latex')
% ylabel('$y$','interpreter','latex')
% set(gca,'Fontsize',18)
% drawnow

