function yprime = pendulum(t,y,FLAG,mu,omegasq);
%  Vdpol(t,y) returns the state derivatives of the Van der Pol equation
%
yprime = [y(2); -omegasq*(y(1)-(mu/6.0)*y(1)^3)];