% 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])