% Marchuk immunology model of E. Hairer, S.P. Norsett, and G. Wanner,
% Solving Ordinary Differential Equations I, Springer-Verlag, Berlin, 
% 1987, p. 297 ff.
%
% There is a state-dependent coefficient xi(y_4(t)) that is not smooth 
% when y_4(t) = 0.1.  The parameter state = +1 if y_4(t) <= 0.1, and -1 
% otherwise.  y_4(t) - 0.1 = 0 is a terminal event.  The integration is 
% restarted with a change of sign in state.  With y_4(0) = 0, state is
% initialized to +1.
%
% There is a discontinuity in the derivative of the history function that
% is handled with 'Jumps'.
%
% There is a parameter h6.  Hairer, Norsett, and Wanner present two plots 
% in Figure 15.8 for h6 = 10 and h6 = 300.  Here the figure for h6 = 300 is 
% reproduced.  This requires scaling of the output and specification of axis.
% The tolerances need to be more stringent than the defaults.

h6 = 300;
state = +1;
opts = ddeset('Jumps',-1e-6,'Events','exer7e','RelTol',1e-5,'AbsTol',1e-8);
sol = dde23('exer7f',0.5,'exer7h',[0, 60],opts,h6,state);
while sol.x(end) < 60
   fprintf('Restart at %4.1f.\n',sol.x(end));
   state = -state;
   sol = dde23('exer7f',0.5,sol,[sol.x(end), 60],opts,h6,state);
end
yplot = sol.y;
yplot(1,:) = 1e4*yplot(1,:);
yplot(2,:) = yplot(2,:)/2;
yplot(4,:) = 10*yplot(4,:);
plot(sol.x,yplot)
legend('1e4*y_1','y_2/2','y_3','10*y_4',0)
title(['Exercise 7.  Marchuk immunology model with h6 = ',num2str(h6),'.'])
xlabel('Restart each time y_4(t) = 0.1.');
axis([0 60 -1 15.5])