function sol = prob_4_8
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

state = +1;

options = ddeset( 'RelTol', 1e-5, 'AbsTol', 1e-8, 'Jumps', [ -10^(-6) ], 'Events', @events );

h6 = 10; 
sol = dde23(@ddes,[0.5 1], @history , [0, 60], options, state, h6 );
while( sol.x < 60 ) 
  if( sol.ie(end) == 1 ) % Integration ended because of an event
    state = -state;
    sol = dde23(@ddes,[0.5 1], sol, [sol.x(end), 60], options, state, h6 );
   end
end
yplot = sol.y;
yplot(1,:) = (1e4)*yplot(1,:); yplot(2,:) = 0.5*yplot(2,:); yplot(4,:) = 10*yplot(4,:); 
figure; plot(sol.x,yplot,'-o'); grid on;  title( [ 'h6 = ',num2str(h6) ] ); 

h6=300;
sol = dde23(@ddes,[0.5 1], @history , [0, 60], options, state, h6 );
while( sol.x < 60 ) 
  if( sol.ie(end) == 1 ) % Integration ended because of an event
    state = -state;
    sol = dde23(@ddes,[0.5 1], sol, [sol.x(end), 60], options, state, h6 );
   end
end
yplot = sol.y;
yplot(1,:) = (1e4)*yplot(1,:); yplot(2,:) = 0.5*yplot(2,:); yplot(4,:) = 10*yplot(4,:); 
figure; plot(sol.x,yplot,'-o'); grid on; title( [ 'h6 = ',num2str(h6) ] ); 


function v = history(t, state, h6 )

v = [ max(0,t+10^(-6)); 1; 1; 0 ]; 


function [value,isterminal,direction] = events(t,y,Z,state,h6)

value = [ y(4)-0.1 ];
isterminal = 1; 
direction = 0; 


function dydt = ddes(t,y,Z,state,h6)

h1 = 2; h2 = 0.8; h3 = 10^4; h4 = 0.17; h5 = 0.5; h7 = 0.12; h8 = 8;
if( state == +1 ) 
  xhi = 1; 
else
  xhi = (10/9)*(1-y(4));
end

dydt = [ (h1-h2*y(3))*y(1);
         xhi*h3*Z(3,1)*Z(1,1) - h5*(y(2)-1);
         h4*(y(2)-y(3)) - h8*y(3)*y(1);
         h6*y(1) - h7*y(4);
         ];