Here I show an alternative code for the solution of the two equations. This code does not use the function ode45 so it might be a little easier for someone to understand how the model works.

%% Initialize variables
b=0.001; % rate of infection
g=0.1; %rate of recovery
i=60; %total run time
S=zeros(i,1); %history of infected over time
I=zeros(i,1); %history of recovered over time
S(1)=499; %population size - 1
I(1)=1; % must start with atleast one infected to spread the disease
time=zeros(i,1); % time history
%% Run model
for t=1:(i-1)
    S(t+1)=-b*S(t)*I(t)+g*I(t)+S(t);
    I(t+1)=b*S(t)*I(t)-g*I(t)+I(t);
    time(t+1)=time(t)+1;
end
%% Display Results
clear g b i t; % delete unnecessary variables
plot(time, I, 'DisplayName', 'I', 'XDataSource', 'time', 'YDataSource', 'I'); hold all; plot(time, S, 'DisplayName', 'S', 'XDataSource', 'time', 'YDataSource', 'S'); hold off; figure(gcf)

An epidemiological graph for the SIS model. S represents the number of susceptibles, I the number of infectious people; When recovered, people become susceptible again. Graph generated using MATLAB code:

P=60
N=500
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
beta = 0.001;
v= 0.1;
[T,Y]=ode45(@sissys1,[0 P],[N 1],options,beta,v);
%plot(T,Y(:,1),T,Y(:,2));
plot(T,Y(:,1),'.',T,Y(:,2),'.');

function dy = sissys1(t,y,beta,v)
dy = zeros(2,1); % a column vector
dy(1) = - beta*y(1)*y(2)+v*y(2);
dy(2) = beta*y(1)*y(2)-v*y(2);
end

I, the copyright holder of this work, release this work into the public domain. This applies worldwide. In some countries this may not be legally possible; if so: I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

Image:Sirsys-p9.png

I, the copyright holder of this work, release this work into the public domain. This applies worldwide. In some countries this may not be legally possible; if so: I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.