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

x=linspace(1.32,2.88,NPT)-2.1;
xdot=linspace(1,2,NPT);
xdotm=linspace(1,2,NPT);
c0=.6;
c1=.66;
for i=1:NPT
    s3=(x(i)-x(i)^3/3-c1)^2/m+k*x(i)^2-c0;
    s2=x(i)-x(i)^3/3-c1;
    s1=m;
    delta=s2^2-s1*s3;
    if(delta>0)
    xdot(i)=(-s2+sqrt(delta))/m;
    xdotm(i)=(-s2-sqrt(delta))/m;
    else
        xdot(i)=-s2;
        xdotm(i)=-s2;
    end
    y(i)=x(i)^2-x(i)^3/3-c1;
end
figure
plot(x,y,'LineWidth',2)
hold on
plot([-1 2],[0,0],'LineWidth',2)
xlabel('x')
ylabel('Fonction x Q(x) ici p=0')
title(' Representation de x Q(x) ici p=0')
%pause
figure
    plot(x,xdot,x,xdotm,'LineWidth',2)% trace de J
   
    xlabel('x')
    ylabel('dx/dt')
     title('Courbe J  definissant un cycle limite')
     %=========