ApproximantCoeffici​entsSIR

This computes the A_n coefficients and Sinf needed for the N-term SIR approximant
given as equation (12a) in Barlow & Weinstein, Physica D 408, 1 (2020) (https://doi.org/10.1016/j.physd.2020.132540)
The inputs correspond to the SIR parameters and initial conditions as
specified by equation (1). Python code available at https://github.com/nsbsma/SIR-approximant

%%% example, reproducing figure 2b in Barlow & Weinstein, Physica D 408, 1 (2020)
alpha=2.73; r=0.0178; I0=7; S0=254; %input parameters
t=0:0.01:6 %time interval 0 to 6 in increments of 0.01
N=25; %number of terms in the approximant, increase until answer stops changing
[A,Sinf,a0,a]=ApproximantCoefficientsSIR(N,alpha,r,S0,I0); %using the code provided to get the stuff needed below
kappa=r*Sinf-alpha;
Ssum=1;
for j=1:N
Ssum=Ssum+A(j)*exp(kappa*t*j); %approximant (equation 4.7a in paper)
end
S=Sinf./Ssum;
plot(t,S,'displayname','S'); hold on
I=alpha/r*log(S/S0)+S0-S+I0;
plot(t,I,'displayname','I');
R=S0+I0-S-I; %b/c S+I+R=constant
plot(t,R,'displayname','R'); xlabel('t'); legend show