function varargout = vdpode(t,y,flag,mu)
%VDPODE Parameterizable van der Pol equation (stiff for large mu).
%   VDPODE(T,Y) or VDPODE(T,Y,[],MU) returns the derivatives vector for the
%   van der Pol equation.  By default, MU is 1, and the problem is not
%   stiff.  Optionally, pass in the MU parameter as an additional parameter
%   to an ODE Suite solver.  The problem becomes more stiff as MU is
%   increased.
%   
%   When MU is 1000 the equation is in relaxation oscillation, and the
%   problem becomes very stiff.  The limit cycle has portions where the
%   solution components change slowly and the problem is quite stiff,
%   alternating with regions of very sharp change where it is not stiff
%   (quasi-discontinuities).  The initial conditions are close to an area of
%   slow change so as to test schemes for the selection of the initial step
%   size.
%   
%   VDPODE(T,Y,'jacobian') or VDPODE(T,Y,'jacobian',MU) returns the Jacobian
%   matrix dF/dY evaluated analytically at (T,Y).  By default, the stiff
%   solvers of the ODE Suite approximate Jacobian matrices numerically.
%   However, if the ODE Solver property Jacobian is set to 'on' with ODESET,
%   a solver calls the ODE file with the flag 'jacobian' to obtain dF/dY.
%   Providing the solvers with an analytic Jacobian is not necessary, but it
%   can improve the reliability and efficiency of integration.
%   
%   VDPODE([],[],'init') returns the default TSPAN, Y0, and OPTIONS values
%   for this problem (see RIGIDODE).  The ODE solver property Vectorized is
%   set to 'on' with ODESET because VDPODE is coded so that calling
%   VDPODE(T,[Y1 Y2 ...]) returns [VDPODE(T,Y1) VDPODE(T,Y2) ...] for
%   scalar time T and vectors Y1,Y2,...  The stiff solvers of the ODE Suite
%   take advantage of this feature when approximating the columns of the
%   Jacobian numerically.
%   
%   L. F. Shampine, Evaluation of a test set for stiff ODE solvers, ACM
%   Trans. Math. Soft., 7 (1981) pp. 409-420.
%   
%   See also ODE15S, ODE23S, ODESET, ODEFILE.

%   Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.10 $  $Date: 1997/11/21 23:27:28 $

if nargin < 4 | isempty(mu)
  mu = 1;
end

switch flag
case ''                                 % Return dy/dt = f(t,y).
  varargout{1} = f(t,y,mu);
case 'init'                             % Return default [tspan,y0,options].
  [varargout{1:3}] = init(mu);
case 'jacobian'                         % Return Jacobian matrix df/dy.
  varargout{1} = jacobian(t,y,mu);
otherwise
  error(['Unknown flag ''' flag '''.']);
end

% --------------------------------------------------------------------------

function dydt = f(t,y,mu)
dydt = [y(2,:); (mu*(1-y(1,:).^2).*y(2,:) - y(1,:))]; % Vectorized

% --------------------------------------------------------------------------

function [tspan,y0,options] = init(mu)
tspan = [0; max(20,3*mu)];              % several periods
y0 = [2; 0];
options = odeset('Vectorized','on');

% --------------------------------------------------------------------------

function dfdy = jacobian(t,y,mu)
dfdy = [ 0                      1
         (-2*mu*y(1)*y(2) - 1)  (mu*(1-y(1)^2)) ];
