a=1;
epsilon=.01;
omega=2;
T=200;
NPT=20000;
dt=T/NPT;
temp=linspace(1,2,NPT)-1;
x=linspace(1,2,NPT);
y=linspace(1,2,NPT);
x(1)=0;y(1)=1;
x(2)=x(1)+dt*y(1);
y(2)=y(1);

for i=3:NPT
    x(i)=2*x(i-1)-x(i-2)-dt*dt*omega*x(i-1)-epsilon*omega*y(i-1)*(x(i-1)^2-a^2)*dt;
    y(i)=(x(i)-x(i-1))/dt;
end
figure
plot(x,y)
title('Un cycle limite de Van der Pol surlequel vient s enrouler la solution')
xlabel('x')
ylabel('dx/dt')