NPT=10000;
T=200;
dt=T/NPT;
m=1;
k=1;

x=linspace(1,3,NPT);
xdot=linspace(1,2,NPT);
xdotm=linspace(1,2,NPT);
for i=1:NPT
    s=k*(.5+2*x(i)-2-x(i)^2/2)/m;
    xdot(i)=sqrt(s);
    xdotm(i)=-xdot(i);
end
figure
    plot(x,xdot,x,xdotm,'LineWidth',2)% trace de H
    hold on
    plot([-1 -1],[-2 2],'-','LineWidth',2)%q(x)=0
    hold on
    plot([1 1],[-2 2],'-','LineWidth',2)% q(x)=0
    hold on
    plot([-2 4],[0 0],'-','LineWidth',2)%axe des x
    hold on
    plot([0 0],[ -2 2],'-','LineWidth',2)
   
    xlabel('x')
    ylabel('dx/dt')
     title('Courbes R- et H ')
     %=========
   % pause
    x(1)=.0;
    x(2)=0;
    xdot(1)=0;
    xdot(2)=0.;
    
    for j=3:NPT
        x(j)=2*x(j-1)-x(j-2)-k*(x(j-1)-2)*dt^2+dt*(1-x(j-1)^2)*(x(j-1)-x(j-2));
        xdot(j)=(x(j)-x(j-1))/dt;
    end
   figure
   plot(x,xdot,'LineWidth',2)
   title('Solution de VDP decalee avec point limite')
   xlabel('x')
   ylabel('dx/dt')