function sol = prob_4_7(r)
% 
% 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.
% 
%-----

c = 1/sqrt(2);
mu = r/10;
options = ddeset( 'RelTol', 1e-5, 'AbsTol', 1e-8, 'Jumps', [ 1-c, 1, 2-c ] );

sol = dde23(@ddes,[1], @history , [0, 10], options, r, mu, c);
figure; plot(sol.x,sol.y(1,:),'-o'); grid on;
title( 'x v.s. y' ); 
figure; plot( sol.x, -sol.yp./(r*sol.y), '-o' ); grid on; axis tight;
title( 'x v.s. yp/y ratio' ); 
tmp = deval( sol, 10 );
fprintf( 'y(10) = %f\n', tmp ); 

function v = history(t, r, mu, c)

v = [ 10 ]; 

function dydt = ddes(t,y,Z,r,mu,c)

if( 0 <= t & t <= 1-c )
  dydt = [ -r*y*0.4*(1-t) ];
elseif( 1-c < t & t <= 1 )
  dydt = [ -r*y*(0.4*(1-t) + 10 - exp(mu)*y) ];
elseif( 1 < t & t <= 2-c )
  dydt = [ -r*y*(10-exp(mu)*y) ];
else
  dydt = [ -r*exp(mu)*y*(Z-y) ];
end