function sol = prob5
% prob5, prob5f of tutorial.

global a b d d1 e gamma T

a  = 1;
b  = 80;
d  = 1;
d1 = 1;
e  = 10;
gamma = 0.5;
T = 0.2;

lambda = 12;
sol = dde23(@prob5f,T,[2; 3.5],[0, 25],[],lambda);

plot(sol.x,sol.y)
legend('I(t)','N(t)',0)
title(['Problem 5. Cooke et al. Epidemic Model with ' ...
       ' \lambda = ',num2str(lambda),'.'])
xlabel('time t')

%=============================================================
function yp = prob5f(t,y,Z,lambda)
global a b d d1 e gamma T
I = y(1);
N = y(2);
Nlag = Z(2,1);
dIdt = lambda*(N - I)*(I/N) - ( d + e + gamma)*I;
dNdt = b*exp(-a*Nlag)*Nlag*exp(-d1*T) - d*N - e*I;
yp = [dIdt; dNdt];