function ex4
%  The Robertson problem is a classic test problem for codes that solve
%  stiff ODEs.  The problem is 
%
%       y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)
%       y(2)' =  0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
%       y(3)' =   3e7*y(2)^2
%
%  It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3) = 0 
%  to steady state.  
%
%  The differential equations satisfy a linear conservation law that can
%  be used to reformulate the problem as the DAE
%
%       y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)
%       y(2)' =  0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
%          0  =  y(1) + y(2) + y(3) - 1
%
%  This problem is used as an example in the prolog to LSODI [1].  Though
%  consistent initial conditions are obvious, the guess y(3) = 1e-3 is
%  used to test initialization.  A logarithmic scale is appropriate for 
%  plotting the solution on the long time interval.  y(2) is small and its 
%  major change takes place in a relatively short time.  Accordingly, the 
%  prolog to LSODI specifies a much smaller absolute error tolerance on
%  this component. Also, when plotting it with the other components, it 
%  is multiplied by 1e4. The natural output of the code does not show 
%  clearly the behavior of this component, so additional output is 
%  specified for this purpose.
%
%  Before effective codes for stiff ODEs were widely available, chemists
%  coped with stiffness by resorting to steady state approximations. The
%  Robertson problems is treated as an example in [2].  The resulting
%  non-stiff DAE is
%
%       y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)
%         0   =  0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
%       y(3)' =   3e7*y(2)^2
%
%  Notice that for y(1) = 1, y(3) = 0, there are two sets of consistent
%  initial conditions.  A guess of y(2) = 1e-3 is used to test initialization.
%  It is found that the steady state assumption works rather well except near
%  the initial point.
% 
%  [1]  A.C. Hindmarsh, LSODE and LSODI, two new initial value ordinary
%       differential equation solvers, SIGNUM Newsletter, 15 (1980), 
%       pp. 10-11.
%
%  [2]  L. Edsberg, Integration package for chemical kinetics, pp. 81-
%       94 in Stiff Differential Systems, R.A. Willoughby, ed., Plenum
%       Press, New York, 1973.

% Copyright 2002, The MathWorks, Inc.

tspan = [0 4*logspace(-6,6)];

% Robertson Problem: Stiff ODE.
problem = 1;
y0 = [1 0 0];  

M = diag([1 1 1]);
% The tolerances are those of the LSODI example.
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);

[t,y] = ode15s(@ex4ode,tspan,y0,options,problem);

y(:,2) = 1e4*y(:,2);
figure(1)
semilogx(t,y)
ylabel('1e4 x y(:,2)')
title('Robertson Problem: Stiff ODE.')

% Robertson Problem: DAE from Conservation Law.
problem = 2;
y0 = [1 0 1e-3];  

M = diag([1 1 0]);
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);

[t,y] = ode15s(@ex4ode,tspan,y0,options,problem);

y(:,2) = 1e4*y(:,2);
figure(2)
semilogx(t,y)
ylabel('1e4 x y(:,2)')
title('Robertson Problem: DAE from Conservation Law.')
xlabel('This is equivalent to the stiff ODEs.')

% Robertson Problem: DAE from Steady State Approximation.
problem = 3;
y0 = [1 1e-3 0];  

M = diag([1 0 1]);
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);

[t,y] = ode15s(@ex4ode,tspan,y0,options,problem);

y(:,2) = 1e4*y(:,2);
figure(3)
semilogx(t,y)
ylabel('1e4 x y(:,2)')
title('Robertson Problem: DAE from Steady State Approximation.')
xlabel('This is only an approximation to the stiff ODEs.')

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

function yp = ex4ode(t,y,problem)

switch(problem)
 case 1

  yp = [ -0.04*y(1) + 1e4*y(2)*y(3)
          0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
          3e7*y(2)^2 ];
  
 case 2
  
  yp = [  -0.04*y(1) + 1e4*y(2)*y(3)
           0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
           y(1) + y(2) + y(3) - 1 ];
  
 case 3
  
  yp = [  -0.04*y(1) + 1e4*y(2)*y(3)
           0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
           3e7*y(2)^2 ];
   
end
