(*******************************************************)
(* Nonlinear oscillator with cubic nonlinear damping
y'' + y = ϵ (y' - (1/3)(((y'))^3)) = 0;
   y(0) = 1;  y'(0) = 0                               *)
(*******************************************************)

Remove["Global`*"]
(* Exact solution *)
ϵ= 0.1;
α=1;
τ=ϵ t;
s=NDSolve[{y''[t] + y[t] == ϵ (y'[t]-(1/3) ((y'[t])^3)),y[0]==0,y'[0]==α},y[t],{t,0,200}]
ypert=  ((α Sin[t])/(Sqrt[(E^(-τ))+((α^2)/4) (1-(E^(-τ)))]));
Dypert = D[ypert,t];

Plot[{Evaluate[y[t]/.s],ypert},{t,0,200},PlotStyle->Thick]
ParametricPlot[{ypert,Dypert},{t,0,200},PlotStyle->Thick]   (* This shows the limit cycle at t->∞ as y=2 Sint, i.e. a circle of radius 2 *)

Remove["Global`*"]
(* Limit cycle from two different initial conditions *)
ϵ= 0.2;
τ=ϵ t;
α1=1;
ypert1=  ((α1 Sin[t])/(Sqrt[(E^(-τ))+(((α1)^2)/4) (1-(E^(-τ)))]));
Dypert1 = D[ypert1,t];

α2 = 3;
ypert2=  ((α2 Sin[t])/(Sqrt[(E^(-τ))+(((α2)^2)/4) (1-(E^(-τ)))]));
Dypert2 = D[ypert2,t];

ParametricPlot[{{ypert1,Dypert1},{ypert2,Dypert2}},{t,0,200},PlotStyle->Thick]