%
% program to solve the single-population epidemic problem
% governing ODE's:
%
%   dS/dt = alpha S - beta IS
%
%   dI/dt = beta IS - gamma I
%
%   dR/dt = gamma I
%
%
%   S(t) the number of susceptibles in the population
%   I(t) the number of infected/infectious individuals in the population
%   R(t) the number of individuals who are immune (through recovery or
%        natural immunity)
%

global alpha beta gamma

figure(1)
clf                           

alpha = 0.02;
beta = 0.001;
gamma = 0.2;


alpha=input('enter alpha (birth-rate): ');
beta=input('enter beta (infection rate): ');
gamma=input('enter gamma (immunity rate): ');



S_0 = 600;
I_0 = 30;
R_0 = 0;

yin(1) = S_0; 
yin(2) = I_0;
yin(3) = R_0;

ratio=S_0*beta/gamma;
fprintf(['epidemiological threshold = ' num2str(ratio) '\n']);



t0 = 0;
tfinal = input(' final time: ');
delt = tfinal/200;
tspan=t0:delt:tfinal;

[t,yode] = ode45('epidemic_rhs',tspan,yin); % variable step RK4/5

S = yode(:,1);
I = yode(:,2);

figure(1)

plot(t,S,'r',t,I,'b')
legend('S(t)','I(t)')

set(gca,'FontName','Helvetica','FontSize',14);
xlabel('time');
ylabel('population');
title(['epidemic: \alpha = ' num2str(alpha), ', \beta = ' num2str(beta), ', \gamma = ' num2str(gamma)]);