% m26chebfun.m - van der Pol oscillator in Chebfun

  C = input('C? (e.g., 1, 10, 0, 0.3) ');
  N = chebop(0,50);
  N.op = @(t,u) diff(u,2) + C*(u^2-1)*diff(u) + u;
  N.lbc = [0.01; 0];
  u = N\0;                                  % solve IVP

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

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

% length(u)
% [val,pos] = max(u,'local')
% diff(pos)