# Epidemic processes with vaccination and immunity loss studied with the BLUES
function method.

Jonas Berx  Institute for Theoretical Physics, KU Leuven, B-3001 Leuven,
Belgium  Joseph O. Indekeu111Given his role as Editor of this journal, Joseph
O. Indekeu had no involvement in the peer-review of articles for which he was
an author and had no access to information regarding their peer-review. Full
responsibility for the peer-review process for this article was delegated to
another Editor. Institute for Theoretical Physics, KU Leuven, B-3001 Leuven,
Belgium

(August 27, 2024)

###### Abstract

The Beyond-Linear-Use-of-Equation-Superposition (BLUES) function method is
extended to coupled nonlinear ordinary differential equations and applied to
the epidemiological SIRS model with vaccination. Accurate analytic
approximations are obtained for the time evolution of the susceptible and
infected population fractions. The results are compared with those obtained
with alternative methods, notably Adomian decomposition, variational iteration
and homotopy perturbation. In contrast with these methods, the BLUES iteration
converges rapidly, globally, and captures the exact asymptotic behavior for
long times. The time of the infection peak is calculated using the BLUES
approximants and the results are compared with numerical solutions, which
indicate that the method is able to generate useful analytic expressions that
coincide with the (numerically) exact ones already for a small number of
iterations.

##  I Introduction

Systems of coupled differential equations (DEs) are omnipresent in the
sciences, medicine and even in the humanities. While some of these systems
have exact closed-form solutions, the majority do not and hence need to be
solved by numerical means. One way to supplement direct numerical simulation
with useful information about the (physical) structure behind the coupled DEs
is to use approximate (semi-)analytical solutions, in which the “physical”
role played by the parameters of the DE is conspicuous.

A stronger need for analytical techniques arises when the parameter space to
be scanned is large. For each set of parameters the numerical computation must
be repeated, which, for large phase spaces, is computationally expensive.
Hence, analytical solutions or approximations are of significant value in this
context. Furthermore, numerical techniques can often be unreliable as a result
of rounding errors, discretization, etc. For specific problems, one usually
requires a specific numerical algorithm that leads to useful results.
(Semi-)analytical methods do not rely on the use of discretization,
perturbation or linearization and can handle most types of nonlinearities
where numerical solvers become very complicated, e.g., integro-differential
equations or fractional differential equations.

The best-known examples of methods to calculate such analytical approximants
are the Adomian decomposition method (ADM) Adomian (1994, 1990), the
variational iteration method (VIM) He (2007) and the homotopy perturbation
method (HPM) He (1999). However, a drawback of these methods is the poor
convergence for large values of the independent variable (e.g., for long
times).

A well-known example of nonlinear coupled DEs is the susceptible-infected-
recovered-susceptible (SIRS) model for studying the time-evolution of
infectious diseases Kermack and McKendrick (1927). There have been attempts to
generate approximate solutions for this model using the ADM Makinde (2007),
the VIM and the HPM Ghotbi  _et al._ (2011); Mungkasi (2021) but these
solutions diverge quickly and necessitate an improvement of the methods in
order to extend the region of convergence.

This paper is structured as follows: in Section II we extend the by now
established BLUES function method Indekeu and Müller-Nedebock (2018); Berx and
Indekeu (2019, 2020, 2021) to a system of nonlinear coupled ordinary DEs.
Next, in Section III, we describe the SIRS model with a constant vaccination
strategy and reduce the system to the two-dimensional subsystem we
subsequently solve. In Section IV we set up the BLUES function method for
coupled DEs and apply it, generating approximate solutions to the SIRS model.
We stress the importance of the identification of the steady states to
construct the associated linear operator. The results are compared with the
numerical solution and with results from three other methods, being the ADM,
VIM and HPM. We use the analytical approximants to calculate the time at which
the infection peak occurs and compare with numerical results. Finally, in
Section V we present the conclusions and an outlook.

##  II BLUES function method for coupled differential equations

Here we extend the BLUES iteration originally developed for ordinary DEs
Indekeu and Müller-Nedebock (2018); Berx and Indekeu (2019, 2020) and partial
DEs Berx and Indekeu (2021) to a system of coupled ordinary DEs. The role of
the inhomogeneous source (or sink) term in the context of the ordinary DE will
now be taken over by a vector of sources (or sinks).

Let us start from an $n$-dimensional system of inhomogeneous nonlinear coupled
ordinary DEs that can be written as a nonlinear operator $\mathcal{N}_{t}$
acting on a vector $\boldsymbol{X}(t)$, with source vector
$\boldsymbol{\chi}(t)$,

$\mathcal{N}_{t}\boldsymbol{X}(t)=\boldsymbol{\chi}(t),\;\forall t>0$ (1)

We now judiciously decompose the nonlinear operator $\mathcal{N}_{t}$ into a
linear operator $\mathcal{L}_{t}$, which contains inter alia a first
derivative in time, and no higher time derivatives, and a residual operator
$\mathcal{R}$, i.e., $\mathcal{R}\equiv\mathcal{L}_{t}-\mathcal{N}_{t}$, which
contains the nonlinear part of $\mathcal{N}_{t}$. We add the subscript $t$ to
the linear and nonlinear operators to emphasize their dependence on time.
Defining suitable initial conditions for $t=0$, i.e.,

$\boldsymbol{X}(0)=\boldsymbol{C}\,,$ (2)

completes the description of the system. In the application we consider in
this work, $\mathcal{R}$ does not depend on $t$. Thus the action of the linear
operator on $\boldsymbol{X}$ results in the following associated linear
coupled system

$\mathcal{L}_{t}\boldsymbol{X}=\boldsymbol{X}_{t}-A\boldsymbol{X}=\boldsymbol{\chi},\;\forall
t>0$ (3)

where the subscript denotes derivative w.r.t. time and the same boundary
conditions are imposed as in (1). The elements of the matrix $A$ are constants
in most of the applications we have in mind. We now propose to rewrite the
system of DEs in an equivalent form by incorporating the initial condition
through multiplication of $\boldsymbol{C}$ with a Dirac delta source
$\delta(t)$ located at $t=0$ and by including this term on the right-hand-side
of the inhomogeneous system, i.e.,

$\mathcal{L}_{t}\boldsymbol{X}=\boldsymbol{X}_{t}-A\boldsymbol{X}=\boldsymbol{\chi}+\boldsymbol{C}\delta\equiv\boldsymbol{\psi},\;\forall
t\geq 0$ (4)

where we have combined the external source $\boldsymbol{\chi}$ and the
“initial condition source” $\boldsymbol{C}\delta$ into the combined source
$\boldsymbol{\psi}$. This formulation amounts to resetting the initial
condition to zero, so that $\boldsymbol{X}(t)=0$ for $t\leq 0^{-}$, followed
by a jump in $\boldsymbol{X}$ implied by integrating the DE over the delta
source, so that $\boldsymbol{X}(0^{+})=\boldsymbol{C}$ and $\boldsymbol{X}(t)$
evolves in a continuous manner for $t>0$. The solution of this linear system
(4) is the following convolution integral Zozulya (2017); Godunov (1997)

$\boldsymbol{X}(t)=(G\ast\boldsymbol{\psi})(t)=G(t)\boldsymbol{C}+\int_{\mathbb{R}}G(t-t^{\prime})\boldsymbol{\chi}(t^{\prime})\mathrm{d}t^{\prime},\;\;\forall
t>0$ (5)

where $G(t)$ is the Green function matrix for the inhomogeneous linear system.
This object can be calculated by finding the matrix exponential
$\exp(At)\equiv\mathbb{1}+At+...$, i.e.,

$G(t)=\mathrm{e}^{At}\Theta(t)\,.$ (6)

with $\Theta(t)$ the Heaviside step function (which we take to be unity for
positive times including $t=0$). This Green function matrix solves the linear
system with a delta function unit matrix source, i.e., it is a solution of the
matrix equation

$G_{t}-AG=\delta(t)\mathbb{1},\;\forall t\geq 0.$ (7)

Adopting the BLUES function strategy, a solution to the nonlinear system (1),
rewritten in the equivalent form

$\mathcal{N}_{t}\boldsymbol{X}=\boldsymbol{\chi}+\boldsymbol{C}\delta\equiv\boldsymbol{\psi},\;\forall
t\geq 0$ (8)

is now proposed in the form of a convolution
$\boldsymbol{X}(t)=(B\ast\boldsymbol{\phi})(t)$, in which the function $B(t)$,
named BLUES function, is taken to be equal to the Green function of the chosen
related linear system, i.e., $B(t)=G(t)$ and the new source
$\boldsymbol{\phi}(t)$ is to be calculated by systematic iteration, using the
given (combined) source $\boldsymbol{\psi}(t)$. This procedure starts from the
following implicit equation, which makes use of the action of the residual
operator,

$\begin{split}\boldsymbol{\phi}&=\boldsymbol{\psi}+(\boldsymbol{\phi}-\boldsymbol{\psi})\\\
&=\boldsymbol{\psi}+\mathcal{R}(B\ast\boldsymbol{\phi}).\end{split}$ (9)

To find the solution to the nonlinear system (1), equation (9) can be iterated
to calculate an approximation for $\boldsymbol{\phi}$ in the form of a
sequence in powers of the residual $\mathcal{R}$. This leads to the sequence
$\boldsymbol{\phi}^{(n)}(t)$, with
$\boldsymbol{\phi}^{(0)}(t)=\boldsymbol{\psi}$. By subsequently taking the
convolution product with $B(t)$, approximate solutions
$\boldsymbol{X}_{\psi}^{(n)}(t)$ to (1) can be obtained. Explicitly, the
iteration for these approximants reads

$\boldsymbol{X}^{(n)}_{\psi}(t)=(B\ast\boldsymbol{\phi}^{(n)})(t)=\boldsymbol{X}^{(0)}_{\psi}(t)+\left(B\ast\mathcal{R}\boldsymbol{X}^{(n-1)}_{\psi}\right)(t)\,,$
(10)

where

$\boldsymbol{X}^{(0)}_{\psi}(t)=(B\ast\boldsymbol{\phi}^{(0)})(t)=(B\ast\boldsymbol{\psi})(t)$
(11)

is the zeroth approximant, which is the convolution product of the linear
problem. We now turn to applying the BLUES iteration procedure to the SIRS
epidemiological model for the spreading of infectious diseases.

##  III The SIRS model with constant vaccination

The SIRS model consists of a group of susceptible $(S)$, infected ($I$) and
recovered/immune ($R$) (human) individuals. The total population $N=S+I+R$ can
grow by virtue of a (constant) birth rate $\pi$ and decay through natural
deaths at a rate $\mu$. At birth, the individuals are vaccinated with
probability $p$ and consequently acquire immunity, effectively adding up to
the group of immune individuals. The remainder of the new births are added to
the susceptible group with probability $1-p$. Through contact with infected
individuals, a person can become infected, following a mass-action law $\beta
SI$ with force of infection $\beta I$. The infected can recover with rate
$\gamma$, acquiring (temporary) immunity and moving to the group of recovered
or immune individuals. Finally, immunity can be lost with rate $\xi$, whereby
recovered individuals move back to the susceptible population. Fig.1
illustrates these different processes in a flow diagram. We assume that all
system parameters and populations are positive, and $0\leq p\leq 1$.

Figure 1: SIRS model with vaccination, immunity loss, births and natural
deaths.

The interactions between the different populations can be described by the
following nonlinear system of DEs, in which the prime denotes derivative
w.r.t. time,

$\displaystyle S^{\prime}(t)$
$\displaystyle=N(t)\pi(1-p)-\beta\frac{S(t)I(t)}{N(t)}-\mu S(t)+\xi R(t)$
(12a) $\displaystyle I^{\prime}(t)$
$\displaystyle=\beta\frac{S(t)I(t)}{N(t)}-(\gamma+\mu)I(t)$ (12b)
$\displaystyle R^{\prime}(t)$ $\displaystyle=N(t)\pi p+\gamma
I(t)-(\mu+\xi)R(t).$ (12c)

The time evolution of the total population $N=S+I+R$ can be found by adding
(12a)-(12c)

$N^{\prime}(t)=(\pi-\mu)N(t)\,,$ (13)

which indicates that the population is not constant. $N(t)$ is a nondecreasing
function of $t$ when the birth rate is higher than or equal to the death rate,
$\pi\geq\mu$. To study the relative importance of the various population
fractions, we scale $S,I$ and $R$ by the total population $N$, i.e.,
$s(t)=S(t)/N(t)$, $i(t)=I(t)/N(t)$ and $r(t)=R(t)/N(t)$. This transforms the
system (12) into the following system for the fractions,

$\displaystyle s^{\prime}(t)$ $\displaystyle=\pi(1-p)-\beta s(t)i(t)-\pi
s(t)+\xi r(t)$ (14a) $\displaystyle i^{\prime}(t)$ $\displaystyle=\beta
s(t)i(t)-(\gamma+\pi)i(t)$ (14b) $\displaystyle r^{\prime}(t)$
$\displaystyle=\pi p+\gamma i(t)-(\pi+\xi)r(t)\,,$ (14c)

where $s(t)+i(t)+r(t)=1$, $\forall t\geq 0$. Note that $\mu$ is eliminated by
this transformation. By using the constraint on the population fractions, we
can eliminate $r(t)$ and study the “two-dimensional” invariant system

$\displaystyle s^{\prime}(t)$ $\displaystyle=\pi(1-p)-\beta
s(t)i(t)-(\pi+\xi)s(t)-\xi i(t)+\xi$ (15a) $\displaystyle i^{\prime}(t)$
$\displaystyle=\beta s(t)i(t)-(\pi+\gamma)i(t)\,.$ (15b)

From a stability analysis performed in appendix A, we deduce that the above
system has two globally stable fixed points, which are known exactly: a
disease-free equilibrium (40a)
$\varepsilon_{0}\equiv(s^{*}_{0},i^{*}_{0})=\left(1-\frac{\pi
p}{\pi+\xi},0\right)$ for which the disease is eradicated and an endemic
equilibrium (40b)
$\varepsilon_{e}\equiv(s^{*}_{e},i^{*}_{e})=\left(\frac{\pi+\gamma}{\beta},\frac{\beta((1-p)\pi+\xi)-(\gamma+\pi)(\xi+\pi)}{\beta(\gamma+\pi+\xi)}\right)$
for which the disease persists and keeps circulating through the population.
The final state of the system is characterized by the _vaccination
reproduction number $R_{V}$_

$R_{V}=\frac{\beta\left((1-p)\pi+\xi\right)}{(\pi+\gamma)(\pi+\xi)}\,,$ (16)

which represents the average number of susceptible individuals which are
infected by one sick individual during their infectious period, while a
vaccination program is in long-time use van den Driessche (2017).

Obviously, the endemic equilibrium $\varepsilon_{e}$ can only exist when
$s^{*}_{e}<1$ and $i^{*}_{e}>0$ which means that the vaccination reproduction
number must satisfy $R_{V}>1$. The disease will be fully eradicated whenever
the disease-free equilibrium is the only possible stable fixed point and the
endemic equilibrium does not exist. This happens when the vaccination
probability $p$ is higher than the critical vaccination threshold $p_{c}$,
which can be inferred from equation (16) by setting $R_{V}=1$,

$p_{c}=\left(\frac{\xi}{\pi}+1\right)\left(1-\frac{\gamma+\pi}{\beta}\right).$
(17)

Note that $p_{c}$ may exceed unity, whereas $p$ cannot. Note that one could
introduce an active constant vaccination strategy which would amount to moving
individuals from $S$ to $R$ with a rate $\omega$. This introduces extra terms
into the nonlinear system (12) and results in transforming the term
$(\pi+\xi)s(t)$ into $(\pi+\xi+\omega)s(t)$ in the $s-$channel of the reduced
subsystem (15). The vaccination reproduction number $R_{V}(\omega)$ in this
situation would then be equal to

$R_{V}(\omega)=\frac{\beta\left((1-p)\pi+\xi\right)}{(\pi+\gamma)(\pi+\xi+\omega)}\,,$
(18)

which amounts to rescaling the $\omega=0$ case, i.e.,

$R_{V}(\omega)=\frac{\pi+\xi}{\pi+\xi+\omega}\,R_{V}(0)\,.$ (19)

It is now clear that $\forall\omega>0$, i.e., active vaccination,
$R_{V}(\omega)<R_{V}(0)$. Hence, active vaccination acts as an additional
mechanism to decrease the reproduction number and possibly to suppress the
existence of an endemic equilibrium. For the sake of simplicity, we will
choose $\omega=0$ in this work.

##  IV BLUES function method for the SIRS model

To find solutions of the system (15), we can write it as a nonlinear matrix
equation, as was demonstrated in section II, i.e.,

$\mathcal{N}_{t}\boldsymbol{X}(t)=\boldsymbol{\psi}(t)$ (20)

with $\boldsymbol{X}(t)$ the vector of solutions,

$\boldsymbol{X}(t)=\begin{pmatrix}s(t)\\\ i(t)\\\ \end{pmatrix}$ (21)

and with source vector $\boldsymbol{\psi}(t)=\boldsymbol{\chi}\
+\boldsymbol{C}\delta(t)$, with $\boldsymbol{\chi}$ the (generally time-
dependent) vector of external sources and $\boldsymbol{C}$ the vector of
initial conditions, i.e.,

$\boldsymbol{\chi}=\begin{pmatrix}\chi_{s}\\\ \chi_{i}\\\
\end{pmatrix}\qquad\;\mbox{and}\;\qquad\boldsymbol{C}=\begin{pmatrix}s_{0}\\\
i_{0}\\\ \end{pmatrix}\equiv\begin{pmatrix}s(0)\\\ i(0)\\\ \end{pmatrix}\,.$
(22)

In this work we will choose a constant vaccination strategy at birth, i.e.,
$\chi_{s}$ and $\chi_{i}$ are time-independent. Note that we have included the
initial conditions in the source $\boldsymbol{\psi}$ by multiplication with a
Dirac point source located at $t=0$, as was explained in Section II.

Now we judiciously tailor the linear operator that is congruous with the
asymptotic equilibrium, by rewriting the nonlinear term in (15) so that
already the linear system captures the stable fixed point exactly. This is
done by including the deviations of the population fractions from their
equilibrium values in the (revised) nonlinear term, as follows,

$\displaystyle s^{\prime}(t)$
$\displaystyle=\pi(1-p)-\beta(s(t)-s^{*})(i(t)-i^{*})-(\pi+\xi-\beta
i^{*})s(t)-(\xi+\beta s^{*})i(t)+\xi+\beta s^{*}i^{*}$ (23a) $\displaystyle
i^{\prime}(t)$ $\displaystyle=\beta(s(t)-s^{*})(i(t)-i^{*})-(\pi+\gamma-\beta
s^{*})i(t)+\beta i^{*}s(t)-\beta s^{*}i^{*}\,,$ (23b)

where $s^{*}$ and $i^{*}$ are the elements of the fixed point vector
$\epsilon=(s^{*},i^{*})$ which represents the equilibrium that is reached.
This equilibrium depends uniquely on the value of $R_{V}$. Note that the
refurbished nonlinear term vanishes at the fixed point and represents the
product of the fluctuations in susceptible and in infected fractions relative
to the equilibrium values. This approach captures the correct asymptotic
behavior for long times provided the linear relaxation times for both $s(t)$
and $i(t)$ exist. This is the case for all $R_{V}\neq 1$. For $R_{V}=1$ this
strategy fails because $i(t)$ is then a marginal variable in the linear system
since $i^{\prime}(t)=0$ at linear level, which precludes an approach to the
fixed point. However, a different choice of linear operator will fix this
problem. We will discuss the special (critical) case $R_{V}=1$ separately. The
novelty of the proposed extension of the BLUES function method to coupled
nonlinear systems lies in the above judicious tailoring of the linear system,
which includes both equilibria by construction.

With the calibration chosen as in (23) we proceed to identify the linear
operator,

$\mathcal{L}_{t}\boldsymbol{X}=\boldsymbol{X}_{t}-A\boldsymbol{X}=\boldsymbol{\chi}+\boldsymbol{C}\delta=\boldsymbol{\psi}$
(24)

where the subscript $t$ on $\boldsymbol{X}$ denotes the time derivative, $A$
is the matrix with elements,

$A=\begin{pmatrix}-(\pi+\xi-\beta i^{*})&-(\xi+\beta s^{*})\\\ \beta
i^{*}&-(\pi+\gamma-\beta s^{*})\end{pmatrix}$ (25)

and $\boldsymbol{\chi}$ is the vector with elements

$\boldsymbol{\chi}=\begin{pmatrix}\pi(1-p)+\xi+\beta s^{*}i^{*}\\\ -\beta
s^{*}i^{*}\\\ \end{pmatrix}\,.$ (26)

The (nonlinear) residual operator $\mathcal{R}$ applied to the solution vector
$\boldsymbol{X}$ then takes the form

$\mathcal{R}\boldsymbol{X}=\begin{pmatrix}-\beta(s(t)-s^{*})(i(t)-i^{*})\\\
\beta(s(t)-s^{*})(i(t)-i^{*})\end{pmatrix}\,.$ (27)

Following the procedure outlined in Section II, we construct an iteration
sequence (10) for the solution vector $\boldsymbol{X}(t)$, i.e.,

$\boldsymbol{X}^{(n)}(t)=(B\ast\boldsymbol{\phi}^{(n)})(t)=\boldsymbol{X}^{(0)}(t)+\left(B\ast\mathcal{R}\boldsymbol{X}^{(n-1)}\right)(t)\,,$
(28)

where $B(t)$ is taken to be the matrix Green function $G(t)$ for the linear
problem defined through (24). This $G(t)$ can be found as the inverse of the
fundamental matrix of the matrix of coefficients $A$ or equivalently as the
matrix exponential of $tA$ multiplied by a step function, i.e.,

$G(t)=\mathrm{e}^{tA}\,\Theta(t)\,.$ (29)

For the disease-free equilibrium ($R_{V}<1$) we obtain,

$G(t)=\begin{pmatrix}\mathrm{e}^{-(\pi+\xi)t}&\frac{\xi+(\pi+\gamma)R_{V}}{(\xi-\gamma)+(\pi+\gamma)R_{V}}\left(\mathrm{e}^{-(\pi+\xi)t}-\mathrm{e}^{-(\pi+\gamma)(1-R_{V})t}\right)\\\
0&\mathrm{e}^{-(\pi+\gamma)(1-R_{V})t}\end{pmatrix}\Theta(t).$ (30)

In this (simple) case the Green function matrix is triangular, so that its
eigenvalues are conspicuous on the main diagonal. These eigenvalues contain
the essential “damping” by virtue of the decaying exponentials, with finite
“linear relaxation times” $\tau_{s}=1/(\pi+\xi)$ and
$\tau_{i}=1/((\pi+\gamma)(1-R_{V}))$. Note that $\tau_{i}$ diverges for
$R_{V}\uparrow 1$. In this limit the relaxation to the disease-free
equilibrium becomes “nonlinear”. The damping (for $R_{V}<1$) ensures that the
long-time asymptotics of the approximants, calculated through convolution, are
well behaved.

The general Green function matrix, appropriate for both disease-free and
endemic equilibria, is more involved and reads,

$G(t)=\frac{\mathrm{e}^{-\frac{Lt}{2}}}{2M}\begin{pmatrix}Z_{+}\,\mathrm{e}^{\frac{Mt}{2}}+Z_{-}\,\mathrm{e}^{-\frac{Mt}{2}}&2\left(\mathrm{e}^{-\frac{Mt}{2}}-\mathrm{e}^{\frac{Mt}{2}}\right)(\beta
s^{*}+\xi)\\\
2\left(\mathrm{e}^{\frac{Mt}{2}}-\mathrm{e}^{-\frac{Mt}{2}}\right)\beta
i^{*}&Z_{+}\,\mathrm{e}^{\frac{-Mt}{2}}+Z_{-}\,\mathrm{e}^{\frac{Mt}{2}}\end{pmatrix}\Theta(t)\,,$
(31)

with

$\begin{split}Z_{\pm}&=M\pm K\\\ K&=\gamma-\xi-\beta(s^{*}+i^{*})\\\
L&=\gamma+\xi+2\pi-\beta(s^{*}-i^{*})\\\
M^{2}&=\gamma^{2}+\left[\beta(s^{*}-i^{*})+\xi\right]^{2}-2\gamma(\gamma-K)\,.\end{split}$
(32)

The zeroth approximant (11) is the convolution of the matrix Green function
with the source vector $\boldsymbol{\psi}(t)$, i.e.,

$\begin{split}\boldsymbol{X}^{(0)}(t)&=(G\ast\boldsymbol{\psi})(t)=\int_{\mathbb{R}}G(t-t^{\prime})\left[\boldsymbol{\chi}+\boldsymbol{C}\delta(t^{\prime})\right]\mathrm{d}t^{\prime}\\\
\end{split}$ (33)

and results in the following expressions for the population fraction of
susceptible and infected individuals, respectively,

$\displaystyle s^{(0)}(t)$ $\displaystyle=\frac{i_{0}(\beta
s^{*}+\xi)}{M}\mathrm{e}^{-Lt/2}\left(\mathrm{e}^{-Mt/2}-\mathrm{e}^{Mt/2}\right)+\frac{s_{0}}{2M}\mathrm{e}^{-Lt/2}\left(Z_{-}\,\mathrm{e}^{-Mt/2}+Z_{+}\,\mathrm{e}^{Mt/2}\right)$
$\displaystyle+\frac{2\beta s^{*}i^{*}(\beta
s^{*}+\xi)}{M}\left(\frac{(\mathrm{e}^{-Lt/2}-\mathrm{e}^{-Mt/2})}{L-M}-\frac{(\mathrm{e}^{-Lt/2}-\mathrm{e}^{Mt/2})}{L+M}\right)$
$\displaystyle+\frac{(\beta
s^{*}i^{*}+\pi(1-p)+\xi)\mathrm{e}^{-Lt/2}}{M}\left(\frac{Z_{-}\,(\mathrm{e}^{Lt/2}-\mathrm{e}^{-Mt/2})}{L+M}\right)$
$\displaystyle+\frac{(\beta
s^{*}i^{*}+\pi(1-p)+\xi)\mathrm{e}^{-Lt/2}}{M}\left(\frac{Z_{+}\,(\mathrm{e}^{Lt/2}-\mathrm{e}^{Mt/2})}{L-M}\right)$
(34) $\displaystyle i^{(0)}(t)$ $\displaystyle=\frac{s_{0}\beta
i^{*}}{M}\mathrm{e}^{-Lt/2}\left(\mathrm{e}^{Mt/2}-\mathrm{e}^{-Mt/2}\right)+\frac{i_{0}}{2M}\mathrm{e}^{-Lt/2}\left(Z_{-}\,\mathrm{e}^{Mt/2}+Z_{+}\,\mathrm{e}^{-Mt/2}\right)$
$\displaystyle-\frac{2\beta i^{*}(\beta
s^{*}i^{*}+(1-p)\pi+\xi)\mathrm{e}^{Mt/2}}{M}\left(\frac{(\mathrm{e}^{-Lt/2}-\mathrm{e}^{-Mt/2})}{L-M}\right)$
$\displaystyle+\frac{2\beta i^{*}(\beta
s^{*}i^{*}+(1-p)\pi+\xi)\mathrm{e}^{Mt/2}}{M}\left(\frac{(\mathrm{e}^{-Lt/2}-\mathrm{e}^{-Mt/2})\mathrm{e}^{-Mt}}{L+M}\right)$
$\displaystyle-\frac{\beta
s^{*}i^{*}}{M}\mathrm{e}^{-Lt/2}\left(\frac{Z_{-}\,(\mathrm{e}^{Lt/2}-\mathrm{e}^{Mt/2})}{L-M}+\frac{Z_{+}\,(\mathrm{e}^{Lt/2}-\mathrm{e}^{-Mt/2})}{L+M}\right)\,.$
(35)

Upon inspection of this and higher approximants (not reported analytically
here) we infer that all BLUES approximants, regardless of the number of
iterations ($n\geq 0$), are qualitatively correct asymptotically, for all
$R_{V}\neq 1$, in that they converge exponentially rapidly towards the exact
fixed point values for long times, in contrast with the other methods which
yield divergences.

We proceed to compare graphically the solution of the SIRS model calculated
with the BLUES method with a precise numerical solution and with approximate
solutions obtained by the ADM, the VIM, or homotopy perturbation method (HPM).
In Table 1 the parameters are shown for three different cases, together with
the values for the vaccination reproduction number $R_{V}$ and critical
vaccination threshold $p_{c}$. Depending on the value of $R_{V}$, we indicate
in the last column of Table 1 the equilibrium attained by the system.

Table 1: Parameters and corresponding equilibria for the 3 studied cases in the SIRS model. The vaccination reproduction number $R_{V}$ and critical vaccination threshold $p_{c}$ are also shown. Note that the latter exceeds unity in Case 2, which is physically equivalent to setting it equal to unity. For all cases $s_{0}=0.8$, $i_{0}=0.2$, $\beta=0.8$, $\gamma=0.03$ and $\pi=0.4$. | $\xi$ | $p$ | $p_{c}$ | $R_{V}$ | Equilibrium
---|---|---|---|---|---
Case 1 | 0.1 | 0.9 | 0.5781 | 0.5209 | $\varepsilon_{0}=(0.28,0)$
Case 2 | 0.5 | 0.9 | 1.0406 $\implies$ 1 | 1.1163 | $\varepsilon_{e}=(0.5375,0.0605)$
Case 3 | 0.1 | 0.5781 | 0.5781 | 1 | $\varepsilon_{0}=(0.5375,0)$

### Case 1: small loss of immunity and high vaccination probability.

As a preliminary remark, we mention that for $\xi=0$ (no loss of immunity) the
SIRS reduces to a SIR model with vaccination, which was treated earlier in
Makinde (2007); Ghotbi  _et al._ (2011) by means of the ADM, HPM and VIM. Here
we consider the SIRS model in which the protection offered by vaccination or
post-disease immunity is lost with a small probability $\xi=0.1$ after some
time. When the vaccination probability $p=0.9$ is higher than the critical
vaccination threshold $p_{c}=0.5781$, the disease will eventually die out and
the system will reach the stable disease-free equilibrium for which
$i\rightarrow 0$ and $s\rightarrow 0.28$. This is shown in Fig. 2. As we
already discussed the BLUES method is accurate and captures the fixed-point
values (40a) of $(s^{*}_{0},i^{*}_{0})$ in the equilibrium exactly. We remark
that the approximants generated by the ADM, VIM and HPM diverge uncontrollably
for longer times while the BLUES approximants converge globally for all $t\geq
0$ and in every iteration.

Figure 2: Comparison between the numerical solution (red line), the fifth-
order ADM (orange dot-dashed line), the fifth-order VIM (blue dotted line),
the fifth-order HPM (green dot-dash-dashed line) and the third BLUES
approximant (black, dashed line) for Case 1 of the SIRS model: $\xi=0.1$ and
$p=0.9$. In this case the disease-free equilibrium is reached for
$t\rightarrow\infty$. Note that the numerical solution and the BLUES
approximant are indistinguishable at this resolution.

We also compare the BLUES approximants for different numbers of iteration and
notice that they converge rapidly towards the numerical solution. This is
shown in Fig. 3.

Figure 3: Comparison and convergence of BLUES approximants after zero (dashed
line), one (dot-dashed line) and two (dotted line) iterations. The numerical
solution is also shown (red line). This figure is for Case 1 of the SIRS
model: $\xi=0.1$ and $p=0.9$. The convergence is global and all approximants
(for all $n$) attain the exact fixed point values for $t\rightarrow\infty$.
The asymptotic behavior for long times is an exponential decay. Note that the
numerical solution and the second BLUES approximant are indistinguishable at
this resolution.

### Case 2: high loss of immunity and high vaccination probability.

As a second example, we consider the case in which immunity is more easily
lost $(\xi=0.5)$ and the population is putting in an effort to vaccinate a
larger number of people $(p=0.9)$. We can deduce from the critical vaccination
probability $p_{c}=1.0406$ in Table 1 that even when all civilians are
vaccinated, immunity is lost so quickly that the population always reaches the
endemic equilibrium and the disease cannot be eradicated. The result of a
comparison between the ADM, VIM, HPM and BLUES methods is shown in Fig. 4. We
also compare the BLUES approximants for different numbers of iteration and
observe that they converge rapidly towards the numerical solution. This is
shown in Fig. 5.

Figure 4: Comparison between the numerical solution (red line), the fifth-
order ADM, VIM, and HPM approximants (respectively, blue dotted, orange dot-
dashed and green dot-dash-dashed lines) and the third BLUES approximant
(black, dashed line) for Case 2 of the SIRS model: $\xi=0.5$ and $p=0.9$. In
this case the endemic equilibrium is reached for $t\rightarrow\infty$. Note
that the numerical solution and the BLUES approximant are indistinguishable at
this resolution.  Figure 5: Comparison and convergence of BLUES approximants
after zero (dashed line), one (dot-dashed line), two (dotted line) and three
(dot-dash-dashed line) iterations. The numerical solution is also shown (red
line). This figure is for Case 2 of the SIRS model: $\xi=0.5$ and $p=0.9$. The
convergence is global and all approximants (for all $n$) attain the exact
fixed point values for $t\rightarrow\infty$ indicated by the horizontal blue
lines. The asymptotic behavior for long times is an exponential decay.

### Case 3: nonlinear relaxation at the dynamical critical point.

In this third case, we study the dynamical criticality at $R_{V}=1$. The
population still reaches the disease-free equilibrium asymptotically, but much
more slowly since in the limit $R_{V}\uparrow 1$ the linear relaxation time
diverges. The “linear” exponential relaxation is replaced by a “nonlinear”
algebraic one, with leading behavior for long times proportional to $1/t$.
This can be inferred exactly from an analysis of the asymptotic behavior of
the system of DEs (23), which at $R_{V}=1$ reduces to,

$\displaystyle s^{\prime}(t)$
$\displaystyle=\pi(1-p)-\beta(s(t)-s^{*})i(t)-(\pi+\xi)s(t)-(\xi+\beta
s^{*})i(t)+\xi$ (36a) $\displaystyle i^{\prime}(t)$
$\displaystyle=\beta(s(t)-s^{*})i(t)\,,$ (36b)

with, in this special case, $s^{*}=s^{*}_{0}=s^{*}_{e}$. Inspection of these
DEs allows us to establish that the leading asymptotic behavior is a $1/t$
power-law decay towards the fixed point,

$\displaystyle s(t)$ $\displaystyle=s^{*}-\frac{1}{\beta\,t}+{\cal O}(t^{-2})$
(37a) $\displaystyle i(t)$
$\displaystyle=\frac{\pi+\xi}{\pi+\xi+\gamma}\,\frac{1}{\beta\,t}+{\cal
O}(t^{-2}),$ (37b)

A successful application of the BLUES function method is possible if we
acknowledge that we must reconsider the decomposition of the problem into a
nonlinear and a linear part. This is necessary in view of the divergence of
the “linear relaxation time” $\tau_{i}=1/((\pi+\gamma)(1-R_{V}))$ and the
concomitant vanishing of one of the eigenvalues of the linear matrix $A$ given
by the triangular form (25), which is appropriate for the disease-free
equilibrium. If we stick to this choice there is not enough “damping” in the
convolution products that govern the iteration procedure (10) and the
approximants diverge, similarly to what routinely happens in the ADM, the VIM
and HPM. However, we can recalibrate the linear operator so that both
eigenvalues are non-zero and nevertheless the correct (disease-free) fixed
point is reached for $t\rightarrow\infty$. This is achieved simply by
retaining the original nonlinear term in (15) without refurbishing it. The
matrix of coefficients $A$ of the linear operator and the residual operator
hence become, respectively,

$A=\begin{pmatrix}-\pi-\xi&-\xi\\\ 0&-\pi-\gamma\\\ \end{pmatrix}$ (38)

and

$\mathcal{R}\boldsymbol{X}=\begin{pmatrix}-\beta s(t)i(t)\\\ \beta
s(t)i(t)\end{pmatrix}\,.$ (39)

This linear system ensures a correct approach to the fixed point for $R_{V}=1$
and is a good starting point for the BLUES iteration in this dynamical
critical point. We recover global convergence to the numerically exact
solution, but at a slower pace than in the noncritical case because the
approximants decay exponentially fast for long times, whereas the exact
solution features an algebraic decay. The characteristic time of the leading
exponentials is constant for all approximants but the amplitude, sign and
polynomial prefactors vary as the iteration number $n$ increases. This is how
the iteration sequence attempts to approximate an algebraic decay in the limit
$n\rightarrow\infty$ (see Fig.7).

Figure 6: Comparison between the numerical solution (red line), the fifth-
order ADM, VIM, and HPM approximants (respectively, blue dotted, orange dot-
dashed and green dot-dash-dashed lines) and the fourth BLUES approximant
(black, dashed line). This figure is for Case 3 of the SIRS model: a critical
point with $\xi=0.1$ and $p=p_{c}=0.5781$, implying $R_{V}=1$. Note that the
BLUES approximant is initially very close to, but later deviates somewhat from
the numerical solution. This is due to the difference in the type of
asymptotic decay, which is conspicuous in Fig.7. Figure 7: Comparison and
convergence of the BLUES approximants after zero (dashed line), two (dot-
dashed line) and four (dotted line) iterations. The numerical solution is also
shown (red line) together with its exact asymptotic behavior (a power-law
decay). This figure is for Case 3 of the SIRS model: a critical point with
$\xi=0.1$ and $p=p_{c}=0.5781$, implying $R_{V}=1$. The convergence is global
and all approximants (for all $n$) attain the exact fixed point values for
$t\rightarrow\infty$. However, the asymptotic behavior of the exact solution
for long times is an algebraic decay, whereas the BLUES approximants feature
an exponential decay.

### Time of the infection peak

We now calculate the (dimensionless) time $t=\hat{t}$ at which the peak of the
infection occurs from the BLUES approximants and compare with the numerically
precise values. At the infection peak, $i^{\prime}(\hat{t})=0$, and hence,
using equation (15b), we deduce that $s(\hat{t})=(\gamma+\pi)/\beta$. So,
instead of trying to solve $i^{\prime}(\hat{t})=0$ directly, we can find
$\hat{t}$ from the susceptible population fraction. The results are shown in
Fig. 8. The BLUES function method accurately captures the infection peak time,
both for the disease-free and the endemic equilibria.

Figure 8: Comparison between the numerically precise (dimensionless) time of
the infection peak $\hat{t}$ in the SIRS model and the values calculated using
the $n$th BLUES approximants ($n=0,1,2,3,4$) for Case 1 (red dashed line and
red circles), Case 2 (black dot-dashed line and black squares) and Case 3
(blue dotted line and blue triangles).

##  V Conclusions

In this paper we have presented a twofold advance. The first consists of
qualitative progress in accuracy and convergence of analytic approximations to
solutions of the SIRS model for epidemic spreading. The second is concerned
with a technical and methodological refinement, being the development of a
matrix BLUES function method for coupled nonlinear ordinary differential
equations.

We have shown that the method can be applied to obtain analytic approximants
for the SIRS model with vital dynamics, constant vaccination strategy and loss
of immunity. We have made a detailed comparison of the iteration procedure
with the Adomian decomposition method, the variational iteration method and
homotopy perturbation method. It is found that all methods succeed in
approximating the (numerically) exact solutions locally and, for the BLUES
function method, also globally. The BLUES method generates approximants with a
higher accuracy, for the same number of iterations performed, and, moreover,
is able to capture the asymptotic behavior of the solutions for long times in
both the disease-free and endemic equilibria. In unpublished work Berx (2021)
we have also studied the SEIRS model, which is an extension of the model
treated in this paper.

A prominent strength of the BLUES function method is that it permits the user
to tailor judiciously the linear part of the operator so that its Green
function contains just enough damping to temper, in each iteration, its
convolution with the emergent source in order to capture the correct
asymptotics of the solution of the nonlinear problem. In other words, for an
optimal choice of the linear part and its associated Green function, the BLUES
iterations converge globally to the numerically exact solution of the problem.
This choice can be viewed as a guiding principle in choosing the associated
linear operator. For a suboptimal choice of the linear part local convergence
can still be achieved, but in that case the advantages over other methods are
less pronounced. An illustration of an optimal and suboptimal choice of linear
operator has been presented in the case of the dynamical critical point of the
SIRS model. In general, when using the BLUES function method to generate
analytical approximants to systems of coupled DEs, one should tailor the
linear system in such a way that it includes all of the existing steady states
and at the same time respects the initial condition by means of a source term,
in order to achieve globally convergent results.

With the ongoing COVID-19 pandemic, it has become clear that mathematical
modeling of epidemic processes is crucial to understanding (and possibly also
predicting) the evolution of these viral outbreaks. None of the
epidemiological models have exact closed-form solutions, with the exception of
the SIS model Hethcote (1989) in which no immunity is possible and individuals
return to the susceptible contingent once they recover from the infection.
Hence, approximate solutions are needed to assess the impact of the different
model parameters without resorting to brute force methods such as numerical
simulation. This paper has aimed at showing that the BLUES method is a good
candidate for obtaining such approximate analytic solutions. This application
of the method to coupled ODEs paves the way for applications to more involved
systems such as coupled PDEs Berx and Indekeu (2021); Horowitz and Kardar
(2019), for instance.

##  Appendix A Stability analysis

The system (15) has two fixed points that can be found by a fixed-point
analysis which reveals a disease-free equilibrium $\varepsilon_{0}$ in which
the disease has died out, and an endemic equilibrium $\varepsilon_{e}$ in
which the infected population density reaches a nonzero asymptotic value,
i.e.,

$\displaystyle\varepsilon_{0}$
$\displaystyle=(s^{*}_{0},i^{*}_{0})=\left(1-\frac{\pi p}{\pi+\xi},0\right)$
(40a) $\displaystyle\varepsilon_{e}$
$\displaystyle=(s^{*}_{e},i^{*}_{e})=\left(\frac{\pi+\gamma}{\beta},\frac{\beta((1-p)\pi+\xi)-(\gamma+\pi)(\xi+\pi)}{\beta(\gamma+\pi+\xi)}\right).$
(40b)

Note that the disease-free equilibrium $\varepsilon_{0}$ is independent of the
average contact rate $\beta$. The endemic equilibrium (40b) can now be
simplified to

$\displaystyle\varepsilon_{e}=(s^{*}_{e},i^{*}_{e})=\left(\left(1-\frac{\pi
p}{\pi+\xi}\right)\frac{1}{R^{({\rm
SIRS})}_{V}},\frac{(1-p)\pi+\xi}{\gamma+\pi+\xi}\left(1-\frac{1}{R^{({\rm
SIRS})}_{V}}\right)\right).$ (41)

These two fixed points are globally stable. This means that for an arbitrary
initial condition in the $si$-plane, the trajectory will converge onto one of
these two fixed points. Which fixed point is reached depends on the system
parameters $\pi,\beta,\gamma,\xi$ and $p$.

The global asymptotic stability for the endemic equilibrium can be proven as
follows. First note that the positive quadrant $\mathbb{R}^{2}_{+}$ of the
$si$-plane is not an invariant set of the system (15), i.e., when $s(t)=0$
then $s^{\prime}(t)<0$ for all values $i(t)>(\pi(1-p)+\xi)/\xi$. This can be
resolved by shifting $(s,i)$ to $(\Sigma,i)$, where

$\Sigma(t)=s(t)+\frac{\xi}{\beta}.$ (42)

Hence, the shifted system becomes

$\displaystyle\Sigma^{\prime}(t)$
$\displaystyle=\pi(1-p)-\beta\Sigma(t)i(t)-(\pi+\xi)\Sigma(t)+\frac{\xi}{\beta}(\xi+\pi)+\xi$
(43a) $\displaystyle i^{\prime}(t)$
$\displaystyle=\beta\Sigma(t)i(t)-(\gamma+\pi+\xi)i(t)\,.$ (43b)

It is now easy to see for $\Sigma(t)=0$, now $\Sigma^{\prime}(t)\geq 0$ for
all values of $i(t)$. By shifting the system, the endemic equilibrium
coordinates change as follows

$(\Sigma^{*}_{e},i^{*}_{e})=\left(\frac{\pi+\gamma+\xi}{\beta},\frac{\pi((1-p)\beta-\gamma-\pi)+\xi(\beta-\gamma-\pi)}{\beta(\pi+\gamma+\xi)}\right).$
(44)

All global properties of the system remain invariant under the shifting of the
coordinates, so we can now try to find a Lyapunov function $V(\Sigma,i)$ to
prove global asymptotic stability of the endemic fixed point (44) of the
shifted system (43) and hence also of the original system (15). We can choose
the following Lyapunov function Korobeinikov and Wake (2002)

$\begin{split}V(\Sigma,i)&=\Sigma^{*}_{e}\left(\frac{\Sigma}{\Sigma^{*}_{e}}-\ln{\frac{\Sigma}{\Sigma^{*}_{e}}}\right)+i^{*}_{e}\left(\frac{i}{i^{*}_{e}}-\ln{\frac{i}{i^{*}_{e}}}\right)\end{split}$
(45)

with time derivative

$\begin{split}V^{\prime}(\Sigma,i)&=\partialderivative{V}{\Sigma}\Sigma^{\prime}(t)+\partialderivative{V}{i}i^{\prime}(t)\\\
&=\left(1-\frac{\Sigma^{*}_{e}}{\Sigma(t)}\right)\Sigma^{\prime}(t)+\left(1-\frac{i^{*}_{e}}{i(t)}\right)i^{\prime}(t).\end{split}$
(46)

Now, from (43a) and (43b) it is clear that for the endemic equilibrium the
following holds

$\begin{split}\beta\Sigma^{*}_{e}i^{*}_{e}&=\pi(1-p)-(\pi+\xi)\Sigma^{*}_{e}+\frac{\xi}{\beta}(\xi+\pi)+\xi\\\
&=(\gamma+\pi+\xi)i^{*}_{e}\,.\end{split}$ (47)

Substituting this property into (46) and simplifying gives after some algebra

$\begin{split}V^{\prime}(\Sigma,i)&=\left(1-\frac{\Sigma^{*}_{e}}{\Sigma(t)}\right)\Sigma^{\prime}(t)+\left(1-\frac{i^{*}_{e}}{i(t)}\right)i^{\prime}(t)\\\
&=-\left(\pi(1-p)+\frac{\xi}{\beta}(\xi+\pi)+\xi\right)\left(\frac{\Sigma(t)}{\Sigma^{*}_{e}}\right)\left(1-\frac{\Sigma^{*}_{e}}{\Sigma(t)}\right)^{2}\\\
&\leq 0\,,\end{split}$ (48)

for all values of $\Sigma,i\geq 0$. This concludes the proof that the endemic
equilibrium is globally asymptotically stable. Proving the global stability of
the disease-free fixed point is now trivial. One can repeat the previous
calculations with the Lyapunov function

$\begin{split}V(\Sigma,i)&=\Sigma^{*}_{0}\left(\frac{\Sigma}{\Sigma^{*}_{0}}-\ln{\frac{\Sigma}{\Sigma^{*}_{0}}}\right)+i\,,\end{split}$
(49)

which is the same as (45) with now $\Sigma^{*}_{0}$ instead of
$\Sigma^{*}_{e}$ and $i^{*}_{e}\rightarrow 0$.

## References

  * Adomian (1994) G. Adomian,  _Solving Frontier Problems of Physics: The Decomposition Method_, 1st ed., Fundamental Theories of Physics, Vol. 60 (Springer Netherlands, 1994).
  * Adomian (1990) G. Adomian, “A review of the decomposition method and some recent results for nonlinear equations,” Math. Comput. Model. 13, 17 – 43 (1990).
  * He (2007) J.-H. He, “Variational iteration method – Some recent results and new interpretations,” J. Comput. Appl. Math. 207, 3 – 17 (2007).
  * He (1999) J.-H. He, “Homotopy perturbation technique,” Comput. Methods Appl. Mech. Eng. 178, 257–262 (1999).
  * Kermack and McKendrick (1927) W. O. Kermack and A. G. McKendrick, “A contribution to the mathematical theory of epidemics,” Proc. R. Soc. Lond. A 115, 700–721 (1927).
  * Makinde (2007) O. Makinde, “Adomian decomposition approach to a SIR epidemic model with constant vaccination strategy,” Appl. Math. Comput. 184, 842 – 848 (2007).
  * Ghotbi  _et al._ (2011) A. R. Ghotbi, A. Barari, M. Omidvar, and G. Domairry, “Application of homotopy perturbation and variational iteration methods to SIR epidemic model,” J. Mech. Med. Biol 11, 149–161 (2011).
  * Mungkasi (2021) S. Mungkasi, “Variational iteration and successive approximation methods for a SIR epidemic model with constant vaccination strategy,” Appl. Math. Model. 90, 1 – 10 (2021).
  * Indekeu and Müller-Nedebock (2018) J. O. Indekeu and K. K. Müller-Nedebock, “BLUES function method in computational physics,” J. Phys. A-Math. Theor. 51, 165201 (2018).
  * Berx and Indekeu (2019) J. Berx and J. O. Indekeu, “Analytic iteration procedure for solitons and traveling wavefronts with sources,” J. Phys. A-Math. Theor. 52, 38LT01 (2019).
  * Berx and Indekeu (2020) J. Berx and J. O. Indekeu, “BLUES iteration applied to nonlinear ordinary differential equations for wave propagation and heat transfer,” J. Phys. A- Math. Theor. 54, 025702 (2020).
  * Berx and Indekeu (2021) J. Berx and J. O. Indekeu, “BLUES function method applied to partial differential equations and analytic approximants for interface growth under shear,” Phys. Rev. Research 3, 033113 (2021).
  * Zozulya (2017) V. V. Zozulya, “Green’s matrices for the BVP for the system of ODE of the first order: applications to the beam theories,” Arch. Appl. Mech 87, 1605–1627 (2017).
  * Godunov (1997) S. K. Godunov,  _Ordinary Differential Equations with Constant Coefficients_ , Translations of mathematical monographs, Vol. 169 (American Mathematical Society, 1997).
  * van den Driessche (2017) P. van den Driessche, “Reproduction numbers of infectious disease models,” Infect. Dis. Model. 2, 288 – 303 (2017).
  * Berx (2021) J. Berx,  _Deposition, diffusion and convection: BLUES approximants and some exact results_ , Ph.D. thesis, KU Leuven (2021).
  * Hethcote (1989) H. W. Hethcote, “Three basic epidemiological models,” in  _Applied Mathematical Ecology_ (Springer Berlin Heidelberg, 1989) pp. 119–144.
  * Horowitz and Kardar (2019) J. M. Horowitz and M. Kardar, “Bacterial range expansions on a growing front: Roughness, fixation, and directed percolation,” Phys. Rev. E 99, 042134 (2019).
  * Korobeinikov and Wake (2002) A. Korobeinikov and G. Wake, “Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models,” Appl. Math. Lett. 15, 955 – 960 (2002).