"Representación del ciclo límite de Van der Pol model, es el sistema con las ecuaciones siguientes";
"Ver las soluciones para μ positivos, nulos y negativos";

Clear[x,y,μ]
μ=1;
eq1:={x'[t]==y[t],y'[t]==μ*(1-x[t]^2)*y[t]-x[t],x[0]==0.5,y[0]==0.5}
eq2:={x'[t]==y[t],y'[t]==μ*(1-x[t]^2)*y[t]-x[t],x[0]==2.1,y[0]==2}

"Isoclinas y campo vectorial";
Isoxx=Solve[0==y,y];
Iso11[x_]=Isoxx[[1,1,2]];
Isoyy=Solve[0==μ*(1-x^2)*y-x,y];
Iso22[x_]=Isoyy[[1,1,2]];isoclinas=Plot[{Iso11[x],Iso22[x]},{x,-3,3},PlotRange->{-3,3},PlotStyle->Thick];
vectores=VectorPlot[{y,μ*(1-x^2)*y-x},{x,-3,3},{y,-3,3}];
Show[{isoclinas,vectores},PlotRange->{{-3,3},{-3,3}},AxesLabel->{"x","y"},AspectRatio->1]

"Resolución";
var:={x,y}
tfinal=200;
trange:={t,0,tfinal}
sol1=NDSolve[eq1,var,trange]
sol2=NDSolve[eq2,var,trange]
graph1=ParametricPlot[Evaluate[{x[t],y[t]}/.sol1],{t,0,tfinal},PlotPoints->500,PlotRange->{{-3,3},{-3,3}},PlotStyle->Thick];
graph2=ParametricPlot[Evaluate[{x[t],y[t]}/.sol2],{t,0,tfinal},PlotPoints->500,PlotRange->{{-3,3},{-3,3}},PlotStyle->Thick];
Show[{graph1,isoclinas,vectores},PlotRange->{{-3,3},{-3,3}},AxesLabel->{"x","y"}]
Show[{graph2,isoclinas,vectores},PlotRange->{{-3,3},{-3,3}},AxesLabel->{"x","y"}]