van der Pol equation

x'' (t)
 + ϵ x' (t) (x (((t))^2) - 1) +  
        x (t) == 0

with the perturbation parameter

ϵ = 1/10

The solution with x[0]=0.5, x'[0]=0

tmpsol = 
      NDSolve[{x''[t]  +  1/10 x'[t] ( x[t]^2 - 1)  +  
              x[t] == 0, x[0] == 0.5, x'[0] == 0}, 
        x[t], {t, 0, 100}]; 
    Plot[ x[t]  /.  tmpsol, {t, 0, 100}, 
      PlotStyle -> {RGBColor[1, 0, 0]}]

Needs["Graphics`Animation`"];

Table[
        ListPlot[{{0, Evaluate[x[t]  /.  tmpsol[[1]]]}}, 
          PlotStyle -> {PointSize[0.04]}, 
          PlotRange -> {{-1, 1}, {-3, 3}}, 
          DisplayFunction -> Identity], {t, 0, 100}];

ShowAnimation[%]

The solution with x[0]=4, x'[0]=0

tmpsol = 
      NDSolve[{x''[t]  +  1/10 x'[t] ( x[t]^2 - 1)  +  
              x[t] == 0, x[0] == 4, x'[0] == 0}, 
        x[t], {t, 0, 100}]; 
    Plot[ x[t]  /.  tmpsol , {t, 0, 100}, 
      PlotStyle -> {RGBColor[0, 0, 1]}]

Table[
        ListPlot[{{0, Evaluate[x[t]  /.  tmpsol[[1]]]}}, 
          PlotStyle -> {PointSize[0.04]}, 
          PlotRange -> {{-1, 1}, {-4, 4}}, 
          DisplayFunction -> Identity], {t, 0, 100}];

ShowAnimation[%]

The solution with x[0]=2, x'[0]=0

tmpsol = 
      NDSolve[{x''[t]  +  1/10 x'[t] ( x[t]^2 - 1)  +  
              x[t] == 0, x[0] == 2, x'[0] == 0}, 
        x[t], {t, 0, 100}]; 
    Plot[ x[t]  /.  tmpsol , {t, 0, 100}, 
      PlotStyle -> {RGBColor[0, 1, 0]}]

The phase space plot for the harmonic oscillator

<< Graphics`PlotField`

PlotVectorField[{y, -x}, {x, -4, 4}, {y, -4, 4}, 
      PlotStyle -> {RGBColor[1, 1, 0]}]

tmpsol1 = 
        NDSolve[{x'[t] == y[t], y'[t] + x[t] == 0, 
            x[0] == 2, y[0] == 0}, {x[t], y[t]}, {t, 0, 
            7}];
tmpsol2 = 
        NDSolve[{x'[t] == y[t], y'[t] + x[t] == 0, 
            x[0] == 0.5, y[0] == 0}, {x[t], y[t]}, {t, 0, 
            7}];
tmpsol3 = 
        NDSolve[{x'[t] == y[t], y'[t] + x[t] == 0, 
            x[0] == 4, y[0] == 0}, {x[t], y[t]}, {t, 0, 
            7}];
ParametricPlot[{{x[t], y[t]} /. tmpsol1,  {x[t], y[t]}  /.  
          tmpsol2,  {x[t], y[t]}  /.  tmpsol3}, {t, 0, 7}, 
      AspectRatio -> 1, 
      PlotStyle -> {RGBColor[1, 1, 0], RGBColor[1, 0, 0], 
          RGBColor[0, 0, 1]}]

The phase space plot for the van der Pol equation

PlotVectorField[{y, -1/10 y ( x^2 - 1)  - x}, {x, -4, 
        4}, {y, -4, 4}]

tmpsol = 
      NDSolve[{y'[t]  +  1/10 y[t] ( x[t]^2 - 1)  +  
              x[t] == 0, x'[t] == y[t], x[0] == 0.5, 
          y[0] == 0}, {x[t], y[t]}, {t, 0, 100}]; 
    ParametricPlot[{x[t], y[t]}  /.  tmpsol ,  {t, 0, 100}, 
      AspectRatio -> 1, PlotRange -> {{-4, 4}, {-4, 4}}, 
      PlotStyle -> {RGBColor[1, 0, 0]}]

tmpsol = 
      NDSolve[{y'[t]  +  1/10 y[t] ( x[t]^2 - 1)  +  
              x[t] == 0, x'[t] == y[t], x[0] == 4, 
          y[0] == 0}, {x[t], y[t]}, {t, 0, 100}]; 
    ParametricPlot[{x[t], y[t]}  /.  tmpsol,  {t, 0, 100}, 
      AspectRatio -> 1, PlotRange -> {{-4, 4}, {-4, 4}}, 
      PlotStyle -> {RGBColor[0, 0, 1]}]

Show[%54, %55, %56]

tmpsol = 
      NDSolve[{y'[t]  +  1/10 y[t] ( x[t]^2 - 1)  +  
              x[t] == 0, x'[t] == y[t], x[0] == 2, 
          y[0] == 0}, {x[t], y[t]}, {t, 0, 100}]; 
    ParametricPlot[{x[t], y[t]}  /.  tmpsol ,  {t, 0, 100}, 
      AspectRatio -> 1, PlotRange -> {{-4, 4}, {-4, 4}}, 
      PlotStyle -> {RGBColor[0, 1, 0]}]

The periodic solution with x[0]=2, x'[0]=0

perturb = 
      ParametricPlot[{2 Cos[t]  +  
            1/10 ( 
                3/4 Sin[t]  - 1/4 Sin[3  t]), 
-2 
              Sin[t] +  1/10 ( 3/4 Cos[t]  - 3/4 Cos[3  t])}, {t, 
          0, 7}, AspectRatio -> 1]

tmpsol1 = 
        NDSolve[{y'[t]  +  1/10 y[t] ( x[t]^2 - 1)  +  
                x[t] == 0, x'[t] == y[t], x[0] == 2, 
            y[0] == 0}, {x[t], y[t]}, {t, 0, 
            7}];
tmpsol2 = 
NDSolve[{y'[t]  +  x[t] == 
            0, x'[t] == y[t], x[0] == 2, 
          y[0] == 0}, {x[t], y[t]}, {t, 0, 7}]; 
    ParametricPlot[{{x[t], y[t]}  /.  tmpsol1 ,  {x[t], y[t]}  /.  
          tmpsol2 , {2 Cos[t]  +  
            1/10 ( 
                3/4 Sin[t]  - 1/4 Sin[3  t]), 
-2 
              Sin[t] +  
            1/10 ( 3/4 Cos[t]  - 3/4 Cos[3  t])}},  {t, 0, 7}, 
      AspectRatio -> 1, 
      PlotStyle -> {RGBColor[0, 1, 0], RGBColor[1, 1, 0], 
          RGBColor[0, 0, 0]}]

ϵ=10

tmpsol = 
        NDSolve[{y'[t]  +  10 y[t] ( x[t]^2 - 1)  +  
                x[t] == 0, x'[t] == y[t], x[0] == 0.5, 
            y[0] == 0}, {x[t], y[t]}, {t, 0, 
            45}];
Plot[ x[t]  /.  tmpsol, {t, 0, 45}, 
        PlotStyle -> {RGBColor[1, 0, 0]}];
ParametricPlot[{x[t], y[t]}  /.  tmpsol ,  {t, 0, 45}, 
      AspectRatio -> 1, PlotRange -> {Automatic, {-15, 15}}, 
      PlotStyle -> {RGBColor[1, 0, 0]}]

Table[
        ListPlot[{0, x[t]}  /.  tmpsol   /.  t -> (t/10),  
          AspectRatio -> 1, 
          PlotRange -> {{-1, 1}, {-3, 3}}, 
          PlotStyle -> {PointSize[0.05], RGBColor[1, 0, 0]}, 
          DisplayFunction -> Identity], {t, 0, 45}];

ShowAnimation[%]