Programming for Dynamical System : SIR Model by Octave

% Function file
% the rhs of the SIR system. Saved as sirrhs.m
function xp = sirrhs(t,x)
global a b;
xp= zeros(2,1);
xp(1) = -a*x(1)*x(2) ;
xp(2) = a*x(1)*x(2)-b*x(2) ;
endfunction

% Call File
% SIR model with constant population. Recovered get immunity.
% Susceptibles, S(t)=x(1)
% Infectives, I(t)=x(2)
% Recovered, R(t)
% dS/dt=-a*S*I, dI/dt=a*S*I-b*I, dR/dt=b*I
% Use that the total population N=S+I+R to solve only first two equations
global a b;
a=.01;
b=.1;
%initial conditions
x0=[99;1];
%Total population
N=100;

% time interval
t0 = 0;
tf =60;
% plot the 2-d vector field with arrows.

figure(1)
clf;
[S,I] = meshgrid(0:10:100);
%xp=sirrhs(1,[S,I]);
Sq = -a*S.*I;
Iq = a*S.*I-b*I;
quiver(S,I,Sq,Iq,1);
%print -depsc2 sir_fig1.eps;

[t x] = ode45('sirrhs',[t0 tf],x0);
%Calculate recovered
R= N-x(:,1)-x(:,2);

figure(2)
clf
plot(t,x(:,1),'b-','LineWidth',2,t,x(:,2),'r-','LineWidth',2,t,R,'g-','LineWidth',2);
ylabel('Population Size','fontsize',15);
xlabel('Time','fontsize',15);
legend('Susceptibles', 'Infectives', 'Recovered');
%print -depsc2 sir_fig2.eps;

figure(3)
clf
quiver(S,I,Sq,Iq,1);
hold on;
plot(x(:,1),x(:,2),'b-','LineWidth',2);
ylabel('Infectives','fontsize',15);
xlabel('Susceptibles','fontsize',15);
%legend('Susceptibles', 'Infectives', 'Recovered');
hold off;
%print -depsc2 sir_fig3.eps;
#######################################
#######################################
For theoretical background and output see: www.ssbtr.net --> Publication ---> Programs
----
Durjoy Majumder, Ph.D.
West Bengal State University 
Secretary, SSBTR