# -*- coding: utf-8 -*-
"""
Created on Wed Jan  1 14:30:33 2015
ok
@author: joseorellana
"""
# import basic modules
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import pylab as pyl
import math
import cmath
import os  


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



for i in range(2,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
	
	
	

fig1=plt.figure()            
plt.plot(x,y)
plt.title('Un cycle limite de Van der Pol surlequel vient s enrouler la solution')
plt.xlabel ('x')
plt.ylabel('dx/dt')

plt.show(fig1)