(* Drawing the phase portrait of the van der Pool oscillator  *)

RungeKuttaStep[f_,y_,y0_,dt_]:=Module[{k1,k2,k3,k4},
  k1=dt N[f /. Thread[y->y0]     ];k2=dt N[f /. Thread[y->y0+k1/2]];
  k3=dt N[f /. Thread[y->y0+k2/2]];k4=dt N[f /. Thread[y->y0+k3  ]];
  y0+(k1+2*k2+2*k3+k4)/6];
EulerStep[f_,y_,y0_,dt_]:=y0+dt N[f /. Thread[y->y0]];
Orbit[Method_,f_,y_,y0_,t_,d_]:=Module[{q=y0,s={y0}},
  Do[q=N[Method[ f,y,q,d]];s=Append[s,q],{Round[t/d]}];s];
Orb[z0_]:={Line[Orbit[RungeKuttaStep, {y,(1-x^2) y - x},{x,y},z0,10,0.1]],
           Line[Orbit[RungeKuttaStep,-{y,(1-x^2) y - x},{x,y},z0,10,0.1]]};
PhasePlane[{{x1_,x2_},{y1_,y2_}},n_]:=
  Show[Graphics[Table[Orb[N[{x1+k*(x2-x1)/n,y1+l*(y2-y1)/n}]],{k,0,n},{l,0,n}]],
  DisplayFunction->$DisplayFunction,PlotRange->{{2x1,2x2},{2y1,2y2}}];
PhasePlane[{{-1.4,1.4},{-1.4,1.4}},2];

