function ex6
%  This is the two phase plug flow 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, specifically the choked flow example.  There are two algebraic 
%  variables.  The variables of HLR correspond to those used here as follows:
%     P <--> y(1), u <--> y(2), y <--> y(3).
%
%  Byrne & Hindmarsh and HLR report a singularity at 1.0958e7, which
%  we also find by going out on minimum step size there.  Integrating
%  just to 1.095e7 is rather faster, which is convenient for testing.

% Copyright 2002, The MathWorks, Inc.

%  Consistent initial conditions from HLR:
y0  = zeros(3,1);
y0(1)  =  1.378e8;
y0(2)  =  2.3665271547403188758740152;
y0(3)  =  4.8146652453167371525906527;

%  Change the algebraic components to test initialization.
y0(2:3) = 0.5 * y0(2:3);

L = 3.2188e7;
%tspan = [0 L];
tspan = [0 1.095e7];

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

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

figure(1)
plot(t,1e-7*y(:,1),t,y(:,3),'--')
axis([0.9e6 1.095e7 0 15])
legend('1e-7*y(:,1)','y(:,3)')
title('Two Phase Plug Flow--Choked Flow.')
xlabel('The first component is multiplied by 1e-7.')

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

function yp = ex6ode(t,y,y0)

  yp = zeros(3,1);
  yp(1) = - y(2)^2;

  R = 45.72;
  rho = 0.814;
  mu = 0.098;
  b = 0.345;
  Qco = 1.7153e6;
  Qa = 3.027e5;

  temp1 = pi*sqrt(0.5*R/rho);
  temp2 = sqrt(0.5*R*rho);
  temp3= log(temp2 * (y(3)/mu) * y(2) - 5);

  yp(2) =   temp1*(R - y(3))^2 * y(2) * (2.5*temp3 + 10.5) ...
           - b*Qco - (y0(1)/y(1)) * Qco * (1 - b);

  yp(3) = 2*temp1*y(2)*((2.5*R*y(3) - 1.25*y(3)^2)* ...
           temp3 + 3*R*y(3) - 2.125*y(3)^2 - 13.6*R*mu/(temp2*y(2))) - Qa;
