<< Graphics`PlotField`
<< Graphics`Graphics`

Q 5.2   (1) analytical  DSolve  ϵ=1.5   p.97

eq = y''[t] - ϵ (1 - y[t]^2) y'[t] +  
            y[t] == 0 /. ϵ ->  1.5

Q 5.2   (2) Numerical  NDSolve   p.97

ans = NDSolve[{eq, y'[0] == 0.1, y[0] == 0}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt = y[t] /. ans

zu11 = 
      Plot[yt, {t, 0, 50}, PlotRange ->  {-4, 4}, 
        AspectRatio ->  0.3, PlotPoints ->  500]

ans = NDSolve[{eq, y'[0] == 10, y[0] == 0}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt1 = y[t] /. ans

Plot[yt1, {t, 0, 50}, PlotRange ->  {-4, 4}, 
      AspectRatio ->  0.3]

ans = NDSolve[{eq, y'[0] == -10, y[0] == 0}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt2 = y[t] /. ans

ans = NDSolve[{eq, y'[0] == 10, y[0] == -3}, 
          y[t], {t, 0, 100}, 
          MaxSteps ->  5000][[1]]
yt3 = y[t] /. ans

Plot[yt3, {t, 0, 50}, PlotRange ->  {-4, 4}, 
      AspectRatio ->  0.3]

ans = NDSolve[{eq, y'[0] == -10, y[0] == 3}, 
          y[t], {t, 0, 100}, 
          MaxSteps ->  5000][[1]]
yt4 = y[t] /. ans

Q 5.2   (3) PlotVectorField   p.97

vt = D[yt, t]; vt1 = D[yt1, t]; vt2 = D[yt2, t]; vt3 = D[yt3, t]; 
    vt4 = D[yt4, t];

zu12 = 
      DisplayTogether[
        ParametricPlot[{yt, vt}, {t, 0, 50}, PlotStyle ->  Hue[0.9]], 
        ParametricPlot[{yt1, vt1}, {t, 0, 50}, PlotStyle ->  Hue[0.7]], 
        ParametricPlot[{yt2, vt2}, {t, 0, 50}, PlotStyle ->  Hue[0.7]], 
        ParametricPlot[{yt3, vt3}, {t, 0, 50}, PlotStyle ->  Hue[0.7]], 
        ParametricPlot[{yt4, vt4}, {t, 0, 50}, PlotStyle ->  Hue[0.7]], 
        PlotVectorField[{v, ϵ (1 - x^2) v - 
                x } /. ϵ ->  1.5, {x, -3, 
            3}, {v, -10, 10}]]

Q 5.2   (1) analytical  DSolve  ϵ=0.2   p.97

eq = y''[t] - ϵ (1 - y[t]^2) y'[t] +  
            y[t] == 0 /. ϵ ->  0.2

Q 5.2   (2) Numerical  NDSolve   p.97

ans = NDSolve[{eq, y'[0] == 0.1, y[0] == 0}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt = y[t] /. ans

zu21 = 
      Plot[yt, {t, 0, 50}, PlotRange ->  {-4, 4}, 
        AspectRatio ->  0.3, PlotPoints ->  500]

ans = NDSolve[{eq, y'[0] == 7, y[0] == 0}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt1 = y[t] /. ans

Plot[yt1, {t, 0, 50}, PlotRange ->  {-4, 4}, 
      AspectRatio ->  0.3]

ans = NDSolve[{eq, y'[0] == -7, y[0] == 0}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt2 = y[t] /. ans

ans = NDSolve[{eq, y'[0] == 5, y[0] == -4}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt3 = y[t] /. ans

Plot[yt3, {t, 0, 50}, PlotRange ->  {-4, 4}, 
      AspectRatio ->  0.3]

ans = NDSolve[{eq, y'[0] == -5, y[0] == 4}, 
          y[t], {t, 0, 50}, 
          MaxSteps ->  5000][[1]]
yt4 = y[t] /. ans

Q 5.2   (3) PlotVectorFieldl   p.97

vt = D[yt, t]; vt1 = D[yt1, t]; vt2 = D[yt2, t]; vt3 = D[yt3, t]; 
    vt4 = D[yt4, t];

zu22 = 
      DisplayTogether[
        ParametricPlot[{yt, vt}, {t, 0, 40}, PlotStyle ->  Hue[0.9]], 
        ParametricPlot[{yt1, vt1}, {t, 0, 30}, PlotStyle ->  Hue[0.7]], 
        ParametricPlot[{yt2, vt2}, {t, 0, 3}, PlotStyle ->  Hue[0.7]], 
        ParametricPlot[{yt3, vt3}, {t, 0, 3}, PlotStyle ->  Hue[0.7]], 
        ParametricPlot[{yt4, vt4}, {t, 0, 3}, PlotStyle ->  Hue[0.7]], 
        PlotVectorField[{v, ϵ (1 - x^2) v - 
                x } /. ϵ -> 0.2, {x, -4, 4}, {v, -7, 
            7}]]

summary

Show[GraphicsArray[{{zu11, zu21}, {zu12, zu22}}]]