% m26_vanderpol.m - van der Pol oscillator

  f = @(t,u,C) [ u(2) ; -u(1)-C*(u(1)^2-1)*u(2) ] ;

  u0 = [0.01;0];
  tspan = [0 50];
  C = input('C? (e.g., 1, 10, 0, 0.3) ');
  [t,u] = ode45(f,tspan,u0,[],C);           % solve IVP

  figure(1)                                 % plot function of t
  plot(t,u(:,1),'k')
  grid on, xlabel t, ylabel('u')
  pause

  figure(2)                                 % plot phase plane
  plot(u(:,1),u(:,2),'r'), axis equal
  grid on, xlabel u, ylabel('u''')