Van der Pol Oscillator

The system we are interested in is the following

(dx/dt)=(1/ϵ)(y-((x^3)/3)+x)
(dy/dt)=-ϵx

Null-Cline Analysis

As an example lets set ϵ=.1.

ϵ=.1;

The null-clines can be found analytically:

Solve[(1/ϵ)(y-((x^3)/3)+x)==0,y]
Solve[-ϵ*x==0,x]

Lets plot the vector field:

plot1=VectorPlot[{(1/ϵ)(y-((x^3)/3)+x),-ϵ*x},{x,-2,2},{y,-1.5,1.5}]

Lets plot the null-clines

plot2=Plot[0.1` (-10.` x+3.333333333333333` (x^3)),{x,-2,2},PlotStyle->{Thick,Red,Dashed}]

plot3=ParametricPlot[{0,y},{y,-2,2},PlotStyle->{Thick,Green,Dashed}]

Now lets combine the plots:

plotVec=Show[plot1,plot2,plot3]

Numerical Analysis

Now, lets numerically solve the equation using NDSolve. We can manually change initial conditions to generate the phase portrait. ParametricPlot is used to plot the trajectory.

s=NDSolve[{x'[t]==1/ϵ*(y[t]-(x[t]^3/3-x[t])),y'[t]==-ϵ*x[t],x[0]==y[0]==1},{x,y},{t,100}];
plottemp=ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,100},PlotStyle->{Thick,Purple}]
plot3=Show[plotVec,plot3,plottemp]

Lets plot how x evolves as a function of time.

plottemp=Plot[Evaluate[x[t]/.s],{t,0,100},PlotStyle->{Thick,Purple}]