function [t,S,I,R] = Program_2_1(beta,gamma,S0,I0,MaxTime)
%
% 
%
% RISK_STRUCTURE( beta, gamma, S0, I0, MaxTime)
%      This is the MATLAB version of program 2.1 from page 19 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the simple SIR epidemic without births or deaths.% 
%

% Sets up default parameters if necessary.
if nargin == 0
   beta=1.4247;
   gamma=0.14286;
   S0=1-1e-6;
   I0=1e-6;
   MaxTime=70;
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 I0<=0 
    error('Initial level of infecteds (%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 MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    
if S0+I0>1
    warning('Initial level of susceptibles+infecteds (%g+%g=%g) is greater than one',S0,I0,S0+I0);
end

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

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

% The main iteration 
[t, pop] = ode45(@Diff_2_1,[0:1:MaxTime],[S I R],[],[beta gamma]);

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

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

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



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

beta=parameter(1); gamma=parameter(2);
S=pop(1); I=pop(2); R=pop(3);

dPop=zeros(3,1);

dPop(1)= -beta*S*I;          %dS/dt = - beta*S*I
dPop(2)= beta*S*I - gamma*I; %dI/dt = beta*S*I - gamma*I
dPop(3)= gamma*I;