ApproximantCoeffici​entsSEIR

%This computes the A_n coefficients, B_n coefficients, lambda_1, and f_infinity needed for the N-term SEIR approximant
%given as equation (15) in the preprint found at https://www.researchgate.net/publication/342211425_Analytic_solution_of_the_SEIR_epidemic_model_via_asymptotic_approximant
%The inputs correspond to the number of terms N (must be even), SEIR parameters, and initial conditions as
%specified in equation (1) in the preprint. The code uses padeapprox.m [1], which implements the algorithm of [2].
%
%[1] https://github.com/chebfun/chebfun/blob/master/padeapprox.m
%[2] P. Gonnet, S. Guettel, and L. N. Trefethen, "ROBUST PADE APPROXIMATION
% VIA SVD", SIAM Rev., 55:101-117, 2013.

% %%% example, reproducing figure 1b in Weinstein et. al.
% alpha=0.466089; beta=0.2; gamma=0.1; I0=0.05; S0=0.88; E0=0.07; R0=0; %input parameters
% t=0:0.1:100 %time interval 0 to 100 in increments of 0.1
% N=18; %number of terms in the approximant, increase until answer stops changing
% [A0,A,B,lambda_1,f_infinity]=ApproximantCoefficientsSEIR(N,alpha,beta,gamma,S0,E0,I0); %using the code provided to get the stuff needed below
% tv=t;clear t; syms t;
% f1=A0; f2=1;
% for j=1:N/2
% f1=f1+A(j)*t.^j;
% f2=f2+B(j)*t.^j;
% end
% FA=f_infinity+exp(lambda_1*t).*(f1./f2);
% dFA=diff(FA,t); t=tv; FA=eval(FA); dFA=eval(dFA);
% SA=exp(FA);
% IA=-1/beta*dFA;
% RA=R0-gamma/beta*(FA-log(S0));
% EA=gamma/beta*(FA-log(S0))-SA+dFA/beta+E0+I0+S0;
% plot(tv,SA,'r','displayname','S');hold on
% plot(tv,EA,'m','displayname','E')
% plot(tv,IA,'c','displayname','I')
% plot(tv,RA,'b','displayname','R'); legend show