% stiff.m  using  ode15s
% y1' = y2    y2' = 1000*(1-y1^2)*y2-y1
function stiff()
  [t,y]=ode15s(@vdp1000, [0 3000], [2 0]);
  plot(t,y(:,1),'-o')
  title('Solution with(\mu = 1000) with ODE15');
  xlabel('Time t');
  ylabel('Solution y');
  legend('y_1','y_2')

  function dydt = vdp1000(t, y)
    dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];
  end % vdp1000
end % stiff

