%
% Cooke, van den Driessche and Zou example (1999)
%
clear all
close all
starttime=clock;
% Parameter values from Fig 3(d)
a = 1;
b = 80;
d1 = 0;
d = 1;
T = [0.2, 1.0, 2.4];
opts = ddeset( 'RelTol', 1.0e-7, 'AbsTol', 1.0e-8 );
% Use a structure for the history, to make computing 
% the residual easier, later.  Take a wide enough interval
% that we cover all three delays used here.
hist = dde23( @(t,y,z) 0, [], 3.5, [-5,0], opts );
% Do three integrations for three different delays
for i=1:3,
   dde = @(t,y,z) b*exp(-a*z-d1*T(i))*z - d*y;
   sol(i) = ddesd(dde, T(i), hist, [0,25], opts );
end;
timetaken=etime(clock,starttime);
% Now plot the solutions
figure
plot( sol(1).x, sol(1).y, 'k-', sol(2).x,sol(2).y, 'k--', ...
      sol(3).x, sol(3).y, 'k-.' )
set(gca,'fontsize',16);
xlabel('t','fontsize',16);
ylabel('y','fontsize',16);
% Now compute and plot the residuals
npts = 10001;
t = linspace(0,25,npts);
r  = zeros(3,npts);
for i=1:3,
    dde = @(t,y,z) b*exp(-a*z-d1*T(i)).*z - d*y;
    [yt,dyt] = deval( sol(i), t );
    % Evaluate the solution at the delayed value of t, too
    zt = deval( sol(i), t-T(i) );
    r(i,:) = dyt - dde(t,yt,zt);
end
figure
r = abs(r);
semilogy( t, r(1,:), 'k-', t, r(2,:), 'k--', t, r(3,:), 'k-.' )
set(gca,'fontsize',16);
xlabel('t','fontsize',16);
ylabel('residual','fontsize',16);
axis([0,25,1E-20,1E5]);
set(gca,'YTick',[1E-20,1E-15,1E-10,1E-5,1,1E5]);
% The code produces the exact solution to a problem within
% 1.0e-5 of the stated prroblem.
   