a=.2;
b=1;
c=.2;
omega=1.5;
T=5000;
NPT=50000;
dt=T/NPT;
x=linspace(1,2,NPT);
y=linspace(1,2,NPT);
r=linspace(1,2,NPT);
phi=linspace(1,2,NPT);
x(1)=.0;
y(1)=.001;
r(1)=sqrt(x(1)^2+y(1)^2);
phi(1)=atan(y(1)/x(1));
for i=2:NPT
    r(i)=r(i-1)+dt*(a*r(i-1)-b*r(i-1)^3);
    phi(i)=phi(i-1)+dt*(omega+c*r(i)^2);
    x(i)=r(i)*cos(phi(i));
    y(i)=-r(i)*sin(phi(i));
end
figure
plot(x,y)
title('Solution de l equation resonante d ordre 3 avec cycle limite')
xlabel ('x')
ylabel('dx/dt')