function ex1
%  This is a singular perturbation problem involving Michaelis-Menten
%  kinetics.  It is studied in Chapter 10 of C. Lin and L. Segal, 
%  Mathematics Applied to Deterministic Problems in the Natural Sciences,
%  SIAM, Philadelphia, 1988.  The problem is
%
%          s0' = -s0 + (s0 + kappa - lambda) * c0    s0(0) = 1   
%  epsilon*c0' =  s0 - (s0 + kappa) * c0             c0(0) = 1
%
%  The parameters used here are taken from p. 317:
%         kappa = 1, lambda = 0.375  and    0 <= t <= 1.2.
%
%  First we write the equations in the form M(t)y' = f(t,y):
%  Let y(1) = s0(t) and y(2) = c0(t).  M(t) = diag{1, epsilon}.  
%
%  This is an ODE for epsilon > 0.  As epsilon -> 0, there is a boundary
%  layer at the origin and for t > 0, the solution tends to an "outer 
%  solution" which satisfies a DAE with epsilon = 0.

% Copyright 2002, The MathWorks, Inc.

kappa = 1;
lambda = 0.375;

tspan = [0 1.2];
y0 = [1; 1];

% First solve the ODE.  Use a "large" value of epsilon so 
% that the boundary layer is visible.

epsilon = 0.01;
M = diag([1 epsilon]);
options = odeset('Mass',M);

[t y] = ode15s(@ex1ode,tspan,y0,options,kappa,lambda);

figure(1)
plot(t,y)
axis([0 1.2 0.4 1.1])
title(['Singularly perturbed ODE with epsilon = ',num2str(epsilon),'.'])

% Now solve as a DAE:
epsilon = 0;
M = diag([1 epsilon]);
options = odeset('Mass',M); 

[tdae ydae] = ode15s(@ex1ode,tspan,y0,options,kappa,lambda);

figure(2)
plot(tdae,ydae)
axis([0 1.2 0.4 1.1])
title('Outer Solution from a DAE.')

fprintf('Given initial conditions are      %6.4f %6.4f.\n',y0(1),y0(2))
fprintf('Consistent initial conditions are %6.4f %6.4f.\n',ydae(1,1),ydae(1,2))

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

function yp = ex1ode(t,y,kappa,lambda)
 
yp = [ (- y(1) + (y(1) + kappa - lambda) * y(2))
         (y(1) - (y(1) + kappa) * y(2)) ];
