%  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 1999, The MathWorks, Inc.

clear all
global y0

%  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];
options = odeset('Mass','M','RelTol',1e-6,'AbsTol',1e-6);

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

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.')


