%
clear all

global a b

dt=0.001;

Nt=8*10^5;

xx=zeros(Nt,1);yy=xx;
T=xx;E=xx;

h=1.0e-8;

xx(1)=1+0.1;yy(1)=0;
xx(1)=h;yy(1)=h;
X=[xx(1);yy(1)];



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

    k1=fdot2d(X,t);
    k2=fdot2d(X+dt/2*k1,t+dt/2);
    k3=fdot2d(X+dt/2*k2,t+dt/2);
    k4=fdot2d(X+dt*k3,t+dt);

    X=X+dt/6*(k1+2*k2+2*k3+k4);
    xx(it)=X(1);
    yy(it)=X(2);
    
    
    
end;
E=(yy.^2)/2-(xx.^2)/2+(xx.^4)/4;

figure(6);
plot(T,E,'r','linewidth',2)
xlabel('t');ylabel('E')


    figure(5)
    plot(xx(1:it),yy(1:it),'r',xx(1),yy(1),'ok',xx(it),yy(it),'+k','linewidth',2)
    drawnow
    
hold on


e=0.1;

xr=roots([1/2,-1,2*e]);




% e=0.1;
% 
% x1=linspace(-10,10,1001);
% y1=sqrt(2*e+x1.^2-x1.^4/2);
% ind=(y1(abs(imag(y1))>1.e-13));
% l1=length(ind);
% y1(ind)=-30*ones(l1,1);
% 
% y2=-sqrt(2*e+x1.^2-x1.^4/2);
% y2(ind)=-30*ones(l1,1);
% figure(500)
% plot(x1,y1,'.k',x1,y2,'.k')
% axis([-2 2 -2 2])