The Van der Pol equation is
$$y''-\mu (1-y^2)y'+y=0, \,\, \mbox{where}\,\, \mu ≥ 0. $$
Use the Runge-Kutta Method to plot the limit cycle and some solutions in the phase plane (the $yy'$-plane) both inside and outside of the limit cycle for a few different values of µ.
The best option for you is to use NDSolve. Since you have not provided specific initial conditions, so I choose randomly.
sol[mu_?NumericQ] := First@NDSolve[{y''[t] - mu*(1 - y[t]^2) y'[t] + y[t] == 0, y[0] == 0,
y'[0] == 1}, y, {t, 0, 10}]
ParametricPlot[Evaluate[{y[t], y'[t]} /. sol[#] & /@ Range[0, 1, 0.2]], {t, 0, 10}]
You can also use ParametricNDSolve,
sol = ParametricNDSolveValue[{y''[t] - mu*(1 - y[t]^2) y'[t] + y[t] == 0, y[0] == 0,
y'[0] == 1}, y, {t, 0, 10}, {mu}];
ParametricPlot[Evaluate[{sol[#][t], D[sol[#][t], t]} & /@ Range[0, 2, 0.1]], {t, 0, 10}]
