function [t,S,I,E,R] = Program_2_6(beta,sigma,gamma,mu,S0,E0,I0,MaxTime)
%
% 
%
% RISK_STRUCTURE( beta, gamma, mu, S0, I0, MaxTime)
%      This is the MATLAB version of program 2.6 from page 41 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the SEIR epidemic with equal births and deaths.
% Note we no-longer explicitly model the recovered class.


% Sets up default parameters if necessary.
if nargin == 0
   beta=520/365.0;
   sigma=1/14.0;
   gamma=1/7.0;
   mu=1/(70*365.0);
   S0=0.1;
   E0=1e-4;
   I0=1e-4;
   MaxTime=60*365;
end

% Checks all the parameters are valid
if S0<=0 
    error('Initial level of susceptibles (%g) is less than or equal to zero',S0);
end

if E0<=0 
    error('Initial level of exposed (%g) is less than or equal to zero',E0);
end

if I0<=0 
    error('Initial level of infectious (%g) is less than or equal to zero',I0);
end

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

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

if mu<0 
    error('Birth / Death rate gamma (%g) is less than zero',mu);
end
    
if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    
if S0+E0+I0>1
    warning('Initial level of susceptibles+infecteds (%g+%g=%g) is greater than one',S0,I0,S0+I0);
end

if beta*sigma<(gamma+mu)*(sigma+mu)
    warning('Basic reproductive ratio (R_0=%g) is less than one',beta*sigma/((gamma+mu)*(sigma+mu)));
end

S=S0; E=E0; I=I0; R=1-S-E-I;

% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_2_6,[0 MaxTime],[S E I],options,[beta sigma gamma mu]);

S=pop(:,1); E=pop(:,2); I=pop(:,3); R=1-S-E-I;

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

subplot(3,1,2) 
h=plot(t,E,'-m',t,I,'-r');
legend(h,'Exposed','Infectious');
xlabel 'Time';
ylabel 'Infected'

subplot(3,1,3)
h=plot(t,R,'-k'); 
xlabel 'Time';
ylabel 'Recovereds'


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

beta=parameter(1); sigma=parameter(2); gamma=parameter(3); mu=parameter(4);
S=pop(1); E=pop(2); I=pop(3);

dPop=zeros(3,1);

dPop(1)= mu -beta*S*I - mu*S;
dPop(2)= beta*S*I - sigma*E - mu*E;
dPop(3)= sigma*E - gamma*I - mu*I;