The Method of Averaging

x' (t)  =  ϵ × (cos^2) 
          t  × \[RawEscape]x (t)

e = 1/100;
sol1 = 
        DSolve[ {x[0] == 1, x'[t] ==   + e x[t] Cos[t]^2}, x, 
          t];

\(\(sol2 = 
        DSolve[\ {x[0] == 1, \(x'\)[t] == \ \ \(+\ e\)\ x[t]\ 1/2}, x, 
          t];\)

\)
Plot[{ x[t]  /.  sol1,   x[t]  /.  sol2,  
        1 +  e Cos[t]^2} ,  {t, 0, 40}]

x' (t)  = -ϵ × (cos^2) 
          t  × \[RawEscape]x (t)

e = 1/10; 
    sol1 = DSolve[ {x[0] == 1, x'[t] ==   - e x[t] Cos[t]^2}, 
        x, t];

sol2 = 
        DSolve[ {x[0] == 1, x'[t] ==   - e x[t] 1/2}, x, 
          t];


Plot[{ x[t]  /.  sol1,  e  Cos[t]^2  /.  sol1,  
        x[t]  /.  sol2} ,  {t, 0, 40}]

Plot[{ x[t]  /.  sol1,  
        e x[t] Cos[t]^2  /.  sol1} ,  {t, 0, 40}]

Plot[{ x[t]  /.  sol1,  
        e x [ t] Cos[t]^2  /.  sol1,  
        e x[2 Pi Floor[t/(2  Pi)]] Cos[t]^2  /.  sol1} ,  {t, 0, 
        40}]

van der Pol eq.

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

with the perturbation parameter

ϵ = 1/10

e = 1/10;

Solution with the method of averaging

xav[t_, C_] := 
        2/(1 + C Exp[-e t])^(1/2) Cos[
            t];
vav[t_, 
        C_] := ((C (E^(-t/10)) Cos[t])/(10 (((1 + C (E^(-t/10))))^(3/2)))) - ((2 Sin[t])/(Sqrt[1 + C (E^(-t/10))]))

av = 
      ParametricPlot[{{xav[t, 0], vav[t, 0]}, {xav[t, 1], 
            vav[t, 1]}, {xav[t, -0.5], vav[t, -0.5]}}, {t, 0, 100}, 
        AspectRatio -> 1]

av = 
      ParametricPlot[{{xav[t, 0], vav[t, 0]}}, {t, 0, 6.28}, 
        AspectRatio -> 1]

Numerical Solution

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}]; 
    amp = 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}]; 
    deamp = ParametricPlot[{x[t], y[t]}  /.  tmpsol,  {t, 0, 100}, 
        AspectRatio -> 1, PlotRange -> {{-4, 4}, {-4, 4}}, 
        PlotStyle -> {RGBColor[0, 0, 1]}]

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}]; 
    per = ParametricPlot[{x[t], y[t]}  /.  tmpsol ,  {t, 0, 100}, 
        AspectRatio -> 1, PlotRange -> {{-4, 4}, {-4, 4}}, 
        PlotStyle -> {RGBColor[0, 1, 0]}]

Show[amp, deamp, per]

Show[per]