function [t,X,Y,Z] = Program_2_3(beta,gamma,nu,mu,rho,X0,Y0,N0,MaxTime)
%
% 
%
  % Program_2_3( beta, gamma, nu, mu, rho, X0, Y0, N0, MaxTime)
%      This is the MATLAB version of program 2.3 from page 35 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the SIR model with a probability of mortality, and
% unequal births and deaths. This code assumes Density-
% Dependent Transmission
%

% Sets up default parameters if necessary.
if nargin == 0
   beta=520/365.0;
   gamma=1/7.0;
   nu=1/(70*365.0);
   mu=1/(70*365.0);
   rho=0.5;
   X0=0.2;
   Y0=1e-4;
   N0=1;
   MaxTime=1e5;
end

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

if Y0<=0 
    error('Initial level of infecteds (%g) is less than or equal to zero',Y0);
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 nu<0 
    error('Birth rate gamma (%g) is less than zero',mu);
end

if mu<0 
    error('Death rate gamma (%g) is less than zero',nu);
end

if rho<0 | rho>1 
    error('Mortality probability rho (%g) is not a valid probability',rho);
end
    
if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    
if X0+Y0>N0
warning('Initial level of susceptibles+infecteds (%g+%g=%g) is greater than population size (%g)',X0,Y0,X0+Y0,N0);
end

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

X=X0; Y=Y0; Z=N0-X-Y;

% The main iteration 
options = odeset('RelTol', 1e-6);
[t, pop]=ode45(@Diff_2_3,[0 MaxTime],[X Y Z],options,[beta gamma nu mu rho]);

X=pop(:,1); Y=pop(:,2); Z=pop(:,3);

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

subplot(3,1,2) 
h=plot(t,Y,'-r');
xlabel 'Time';
ylabel 'Infectious'

subplot(3,1,3)
h=plot(t,Z,'-k', t,X+Y+Z,'--k'); 
xlabel 'Time';
ylabel 'Recovereds (-), Total Population (--)'


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

beta=parameter(1); gamma=parameter(2); nu=parameter(3); mu=parameter(4); rho=parameter(5);
X=pop(1); Y=pop(2); Z=pop(3);

dPop=zeros(3,1);

dPop(1)= nu -beta*X*Y - mu*X;
dPop(2)= beta*X*Y - (gamma+mu)*Y/(1-rho);
dPop(3)= gamma*Y - mu*Z;