function [t,S,I,R] = Program_5_4(beta0, beta1 ,gamma,alpha0, alpha1, mu, m, K, X0,Y0,N0,MaxTime)
%
% 
%
% Program_5_4( beta, gamma, alpha0, alpha1, mu, m, K, X0, Y0, N0, MaxTime)
%      This is the MATLAB version of program 5.1 from page 186 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is a model for Rabit Hemorrhagic Disease, in which both transmission
% rate and birth rates can be seasonally forced.
% Note that we are using numbers of individuals in this formulation.
% Bifurcation plots are not possible with this code.
%

% Sets up default parameters if necessary.
if nargin == 0
   beta0=0.936;
   beta1=0.1;
   gamma=0.025;
   alpha0=0.02;
   alpha1=0.1;
   mu=0.01;
   m=0.475;
   K=10000;
   
   X0=0.5;
   Y0=0.01;
   N0=0.6;
   MaxTime=60*365;
end

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

if Y0<=0 
    error('Initial number of infecteds (%g) is less than or equal to zero',Y0);
end

if N0<=0 
    error('Initial population size (%g) is less than or equal to zero',N0);
end

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

if beta1<0 
    error('Seasonality beta1 (%g) is less than zero',beta1);
end

if beta1>=1
    error('Seasonality beta1 (%g) is greater than or equal to one',beta1);
end
    
if alpha1<0 
    error('Seasonality alpha1 (%g) is less than zero',alpha1);
end

if alpha1>=1
    error('Seasonality alpha1 (%g) is greater than or equal to one',alpha1);
end
    
if gamma<=0 
    error('Recovery rate gamma (%g) is less than or equal to zero',gamma);
end

if alpha0<0 
    error('Birth rate alpha_0 (%g) is less than zero',alpha0);
end

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

if m<0 
    error('Mortality rate due to infection m (%g) is less than zero',m);
end

if K<=0 
    error('Carrying capacity K (%g) is less than or equal to zero',K);
end
    
if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    
if X0+Y0>N0
    error('Initial level of susceptibles+infecteds (%g+%g=%g) is greater than the population size %g',X0,Y0,X0+Y0,N0);
end


% The main iteration 

    options = odeset('RelTol', 1e-5);
    [t, pop]=ode45(@Diff_5_4,[0 MaxTime],[X0 Y0 N0],options,[beta0 beta1 gamma alpha0 alpha1 mu m K]);

    T=t/365; X=pop(:,1); Y=pop(:,2); N=pop(:,3);

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

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

    subplot(3,1,3)
    h=plot(T,N,'-k');
    xlabel 'Time (years)';
    ylabel 'Population Size'



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

beta0=parameter(1); beta1=parameter(2); gamma=parameter(3); 
alpha0=parameter(4); alpha1=parameter(5);
mu=parameter(6); m=parameter(7); K=parameter(8);

alpha=alpha0*(1+alpha1*sin(2*pi*t/365));
beta=beta0*(1+beta1*sin(2*pi*t/365));

X=pop(1); Y=pop(2); N=pop(3);

dPop=zeros(3,1);

dPop(1)= alpha*N -beta*X*Y - (mu+N/K)*X;
dPop(2)= beta*X*Y - (gamma + mu + m + N/K)*Y;
dPop(3)= (alpha - mu - N/K)*N - m*Y;



