%  This is the DAS 2 test problem of I. Cameron, Solution of differential-
%  algebraic systems using diagonally implicit Runge-Kutta methods, IMA
%  J. Num. Anal., 3 (1983) 273-289.  The initial conditions specified by
%  Cameron are not consistent.  There are a number of events leading to
%  discontinous changes in the equations.  These must be solved as fresh
%  problems because the solution must change discontinuously at an event
%  to be consistent.  The effects of some of the events are very much
%  faster than the scale of the problem.  To illustrate this, the action
%  of a valve is displayed after the first event for a time interval 
%  of ten times the time constant of the valve.  The same data is then 
%  displayed on the time scale of the whole simulation to show that it 
%  appears as a discontinuity on this scale.  The final figure for the 
%  height of the fluid shows what appear to be discontinuities at the times
%  of the events, which are indicated on the plot.
%
%  Cameron's variables correspond to those used here as follows:
%        y(1) = y, y(2) = s, y(3) = z, y(4) = T_G, y(5) = A_v,
%        y(6) = P_G, y(7) = F_1, y(8) = F_2, y(9) = V_G, y(10) = P2.

% Copyright 1999, The MathWorks, Inc.

clear all
global I P1 P3
I = 1;
P1 = 400;
P3 = 100;

tau = 2.77e-4;        %  Time constant of valve.

y0 = zeros(10,1);
y0(2) = 1;
y0(4) = 288;
y0(5) = 1;
y0(6) = 100;
y0(7) = 59.154;
y0(9) = 201.0619;
y0(10) = 100;

options = odeset('Mass','M');

tspan = [0 1];
events = 0;
[t,y] = ode15s('ex10ode',tspan,y0,options);

%  Now a change in the value of I:
I = 0.7;
tspan = [1 1.5];
events = [events 1];
nout = length(t);
y0 = y(nout,:);
htry = t(end) - t(end-1);
[tnext,ynext] = ode15s('ex10ode',tspan,y0,odeset(options,'InitialStep',htry));
t = [t; tnext];
y = [y; ynext];
%  Plot the transient in y = y(1):
i = find(tnext < 1+10*tau);
figure
plot(tnext(i),ynext(i,1),'*')
legend('Action of valve ')
xlabel('Time interval of 10*tau after first event.')
figure
plot(tnext(i),ynext(i,1))
axis([0 10 -500 100])
legend('Action of valve ')
xlabel('Same as previous figure on the time scale of the problem.')

%  Now another change in the value of I:
I = 0.6;
tspan = [1.5 2];
events = [events 1.5];
nout = length(t);
y0 = y(nout,:);
htry = t(end) - t(end-1);
[tnext,ynext] = ode15s('ex10ode',tspan,y0,odeset(options,'InitialStep',htry));
t = [t; tnext];
y = [y; ynext];

%  Now a change in the value of P1:
P1 = 500;
tspan = [2 2.5];
events = [events 2];
nout = length(t);
y0 = y(nout,:);
htry = t(end) - t(end-1);
[tnext,ynext] = ode15s('ex10ode',tspan,y0,odeset(options,'InitialStep',htry));
t = [t; tnext];
y = [y; ynext];

%  Now a change in the value of P3:
P3 = 110;
tspan = [2.5 3];
events = [events 2.5];
nout = length(t);
y0 = y(nout,:);
htry = t(end) - t(end-1);
[tnext,ynext] = ode15s('ex10ode',tspan,y0,odeset(options,'InitialStep',htry));
t = [t; tnext];
y = [y; ynext];

%  Now another change in the value of I:
I = 0.7;
tspan = [3 3.5];
events = [events 3];
nout = length(t);
y0 = y(nout,:);
htry = t(end) - t(end-1);
[tnext,ynext] = ode15s('ex10ode',tspan,y0,odeset(options,'InitialStep',htry));
t = [t; tnext];
y = [y; ynext];

%  Now changes in both P1 and P3:
P1 = 400;
P3 = 100;
tspan = [3.5 10];
events = [events 3.5];
nout = length(t);
y0 = y(nout,:);
htry = t(end) - t(end-1);
[tnext,ynext] = ode15s('ex10ode',tspan,y0,odeset(options,'InitialStep',htry));
t = [t; tnext];
y = [y; ynext];

%  All done!
figure
plot(t,y(:,3),'-o')
set(gca,'YLim',[0 1.3e-6]);

figure
plot(t,y(:,3),'-o')
set(gca,'XLim',[1 1.006],'YLim',[2.48e-7 2.64e-7]);


