SIR Simulation: Final Epidemic Size

%% SIR Simulation: Final Epidemic Size
% This code piggy-backs off a code with a similar name but this coed
% can estimate the final epidemic size

clc, clear
Beta = 0.2:0.01:1;
mu = 0.1;
S = NaN(size(Beta));
I = NaN(size(Beta));
R = NaN(size(Beta));
BetaT = NaN(size(Beta));
dt = 1;
t = 0:1:600;

BetaT(1:end) = Beta;
N = 8000000;
R(1) = 0;
I(1) = 1;
S(1) = N-I(1);
stop_length = length(Beta);

for i = 1:stop_length
    S(i+1) = ((-BetaT(i)*I(i)*S(i))/N)*dt+S(i);
    I(i+1) = (((BetaT(i)*I(i)*S(i))/N)-(mu*I(i)))*dt+I(i);
    R(i+1) = (mu*I(i))*dt+R(i);
    Rss1(i) = R(i);
    R0(i) = Beta(i)/mu;
    Rss2(i) = N*(1-exp(-R0(i)));
end


hold on
plot(Beta, Rss1)
plot(Beta, Rss2)
legend('Numerical Solutions for Rss', 'Analytical Solutions for Rss', 'Location', 'best')
ylabel('Epidemic Size')
xlabel('Beta')
title('Numerical and Analytical Solutions for Rss vs. beta')