function sol = prob_4_12
% 
% 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.
% 
%-----

opts = ddeset( 'RelTol', 1e-5, 'AbsTol', 1e-8 );

Delays = [ 0.2 1.0 2.4 ];
B = [ 20 80 20 80 ];
D1 = [ 1 1 0 0 ];
for dataset=1:4,
  a=1;
  d=1;
  d1=D1(dataset);
  b=B(dataset);
  for i=1:3
    T = Delays(i);
    sol(i) = dde23(@ddes,T,3.5,[0,25],opts,a,d,d1,b,T);
   end
   figure; plot( sol(1).x, sol(1).y, '-og', ...
                 sol(2).x, sol(2).y, '-ob', ...
                 sol(3).x, sol(3).y, '-or' ); grid on;
end

return;

function dydt = ddes(t,y,Z,a,d,d1,b,T)

dydt = [ b*Z(1,1)*(exp(-a*Z(1,1)))*exp(-d1*T) - d*y(1);
         ];