function [t,S,I] = Program_3_5(n, m, beta, gamma, mu, S0, I0, MaxTime)
%
% 
%
% Program_3_5( n, m, beta, gamma, mu, S0, I0, MaxTime )
%      This is the MATLAB version of program 3.5 from page 94 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the S(E)IR model with multiple stages to create gamma-distributed
% exposed and infectious periods
%
% n is the total number of infection classes; classes 1 to m are exposed,
% classes m+1 to n are infectious.
% Birth and death rates are assumed equal.
%
% As well as the default, you may want to compare the structured model:
% Program_3_5(10, 0, 1, 0.1, 0, 0.5, 1e-6, 60);
% with the unstructured version
% Program_3_5(1, 0, 1, 0.1, 0, 0.5, 1e-6, 60);
% Or compare the SEIR:
% Program_3_5(10, 5, 1, 0.1, 0, 0.5, 1e-4, 150);
% with the unstructured version
% Program_3_5(2, 1, 1, 0.1, 0, 0.5, 1e-4, 150);


% Sets up default parameters if necessary.
if nargin == 0
   n=13;
   m=8;
   gamma=1/13;
   beta=17/5;
   mu=1/(55*365);
   S0=0.05;
   I0=0.00001;
   MaxTime=30*365;
end

% Checks all the parameters are valid


if I0<=0
    error('Initial level of infection (%g) is less than or equal to zero',I0);
end
if S0<0
    error('Initial level of susceptibles (%g) is less than zero',S0);
end
if n<=0
    error('Number of infected stages (%d) is less than or equal to zero',n);
end
if m>=n
    error('Number of infectious stages (%d) is less than or equal to zero\nm is probably too big',m-n);
end

if gamma<=0
    error('Recovery rate gamma (%g) is less than or equal to zero',gamma);
end

if beta<=0
    error('Transmission rate beta (%g) is less than or equal to zero',beta);
end
      
if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    
I0=I0*ones(1,n)/n;

% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_3_5,[0 MaxTime],[S0 I0],options,[n m beta gamma mu]);

% plots the graphs with scaled colours
subplot(3,1,1)
h=plot(t,pop(:,1),'-g');
xlabel 'Time';
ylabel 'Susceptibles, S'

subplot(3,1,2)
h=semilogy(t,pop(:,2:end),'-r');
set(h(1:m),'Color',[1 0.7 0]);
set(h((m+1):n),'LineWidth',1);
if m>0
    legend(h([1 n]),'Exposed','Infectious');
end
xlabel 'Time';
ylabel 'Infection, I_i'

subplot(3,1,3) 
if n>1
    h=plot(t,sum(pop(:,2:end)'),'-r'); 
else
    h=plot(t,pop(:,2),'-r'); 
end
xlabel 'Time';
ylabel 'Total Infection'



% Calculates the differential rates used in the integration.
function dPop=Diff_3_5(t,pop, parameter)

n=parameter(1); m=parameter(2); beta=parameter(3); gamma=parameter(4); mu=parameter(5);

S=pop(1); I=pop(2:end); Inf=sum(I(m+1:n));

dPop=zeros(n+1,1);

dPop(1) = mu - beta*Inf*S - mu*S;
dPop(2) = beta*Inf*S - gamma*n*I(1) - mu*I(1);
for i=2:n
    dPop(i+1) = gamma*n*I(i-1) - gamma*n*I(i) - mu*I(i);
end