function ex13
%  This is the discharge pressure control problem from E. Hairer, Ch. Lubich,
%  and M. Roche, The Numerical Solution of Differential-Algebraic Systems by
%  Runge-Kutta Methods, Lecture Notes in Math. #1409, Springer, Berlin, 1989.
%  As posed by Preston, Berzins, Dew, and Scales, HLR show that the problem
%  is of index 2.  They do this by reducing it to a problem that is shown
%  to be of index 1.  The index 1 problem is solved here.  The variables of
%  their analysis correspond to the ones used here as follows:
%
%    k = y(1), c = y(2), M = y(3), p = y(4), P = y(5), m = y(6).
%
%  For analytical purposes HLR introduce a variable s = c + p/15.  This is
%  not used here with the consequence that the mass matrix is not diagonal,
%  hence the problem is not in semi-explicit form.  HLR solve the problem
%  with RADAU5, which requires consistent initial conditions.  The values
%  given the code here are inconsistent values that approximate well their
%  consistent values.  The consistent initial conditions computed by the
%  code are close to the values of HLR, but they are not identical.  The
%  plot of the given flow in and the computed flow out computed here agrees 
%  with Fig. 9.8 of HLR.

% Copyright 2002, The MathWorks, Inc.

mu = 15 + 5*tanh(0 - 10);
y0 = zeros(6,1);

%  Use inconsistent initial values approximating well the consistent 
%  initial values used by HLR.
y0(1) = 0.25;
y0(2) = 0.25;
y0(3) = 734;
y0(4) = 99.1;
y0(5) = 36.7;
y0(6) = 9.998;

M = diag([1 1 1 0 0 0]);
M(2,4) = 1/15;

options = odeset('Mass',M);

[t,y] = ode15s(@ex13ode,[0 40],y0,options);

mu = 15 + 5*tanh(t - 10);
plot(t,mu,t,y(:,6))
legend('given flow in ','computed flow out ')
axis([0 40 0 40])
title('Discharge Pressure Control Problem')

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

function yp = ex13ode(t,y)

  %  k = y(1), c = y(2), M = y(3), p = y(4), P = y(5), m = y(6)

  mu = 15 + 5*tanh(t - 10);
  mup = 5/cosh(t - 10)^2;

  f = zeros(6,1);
  f(1) =   (y(2) - y(1))/20;
  f(2) = - (y(4) - 99.1)/75;
  f(3) =    mu - y(6);
  f(4) =   (3.35 - 0.075*y(6) + 0.001*y(6)^2) - y(4)/y(5);
  f(5) =    20*y(5) - y(3);
  f(6) =    y(3)*f(3)/400 + mu*mup/(1.2*y(1))^2 ...
            - mu^2 * f(1)/(1.2^2 * y(1)^3);
  yp = f;
