function [t,y] = sirdelaymodel
% TW Carr for SMU math 6350
% SIR model with temporary immunity

global b mu btrans grecov v tau

% Set the parameter values all in units of 1/yrs
b =  1/75;    			% Birth rate
mu = 1/75;			% Death rate
btrans = ?.?;			% Transmission rate (beta is a predefined matlab function)
grecov = ?.?;			% Recovery rate (gamma is a predefined matlab functin)
% Determine the basic reproductive number
R0 = btrans/(grecov+mu)

% Total integration time
tend = ?;
% Specify the history function (initial conditions)
y0history = [0.99, 0.01, 0.0];	% [S, I, R]
% Specify the delay time
tau = 1.;

% Simulate the DDEs
options = ddeset('RelTol',1e-?,'AbsTol',1e-?);
sol = dde23(@f,[tau],y0history,[0,tend],options);

% Find the last time for which computations were successful
tsize = size(sol.x);
tcomp = sol.x(tsize(2))

% Evaluate the solution for analysis and plotting. ' = transpose.
t = linspace(0,tcomp,1000); 	% Generate n equally space points between 0 and tend.
y = deval(sol,t)';

figure(1)
plot(t,y(:,1),'r',t,y(:,2),'g-',t,y(:,3),'b:');

end

%------------------------------------------------------------------------
function dydt = f(t,y,ydelay)

global b mu btrans grecov tau

dydt = zeros(3,1);

S = y(1);		% Dependent variables
I = y(2);
R = y(3);
SIRdelay = ydelay(:,1);	% Delay varables

dydt(1) = mu*(1.0 - S) - btrans*(S*I) + grecov*exp(-mu*tau)*SIRdelay(2);
dydt(2) = btrans*(S*I) - (mu+grecov)*I;
dydt(3) = grecov*I - mu*R - grecov*exp(-mu*tau)*SIRdelay(2);

end


