tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
[t,y] = ode45(@vanderpoldemo, tspan, y0,[],Mu);

% Plot of the solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')
disp(' \rho')