
function [S,I,R] = disease(N,m,b,g)

% m birth and death of Susceptibles.
% b SPREAD OF DISEASE
% g recovery rate of Infecteds
% N is the number of time steps required


S = zeros(1,N+1);
I = zeros(1,N+1);
R = zeros(1,N+1);

S(1) = 0.9;
I(1) = 0.1;
R(1) = 0;
t = 0;
i=0;

for i=1:N
 S(i+1) = S(i) + m - b*S(i)*I(i) - m*S(i);
 I(i+1) = I(i) + b*S(i)*I(i) - g*I(i) - m*I(i);
 R(i+1) = R(i) + g*I(i) - m*R(i);
end
t = [0:1:N];
if N>50
    plot(t,S,t,I,t,R,'LineWidth',2)
else
    plot(t,S,'kd-',t,I,'ro-',t,R,'bx-','LineWidth',2)
end
legend('Susceptibles, S','Infected, I','Recovered, R','Location','East');
axis([0 N 0 1]);
xlabel('i','FontSize',16);
ylabel('percentage population','FontSize',16);
tit = sprintf('m = %.3f, b = %.3f, g = %.3f',m,b,g);
title(tit,'FontSize',14)
