Accurate Numerical Treatment on a Stochastic SIR Epidemic Model with Optimal Control Strategy

Abstract

In this paper, a numerical study has been undertaken on the susceptible-infected-recovered (SIR) epidemic model that encompasses the mechanisms of the evolution of disease transmission; a prophylactic vaccination strategy in the susceptible populations, depending on the infective individuals. We furnish numerical and graphical simulation combined with explicit series solutions of the proposed model using the New Iterative Method (NIM) and Modified New Iterative Method (MNIM). The analytic-numeric New Iterative Method failed to deliver accurate solution for the large time domain. A new reliable algorithm based on NIM, the coupling of the Laplace transforms, and the New Iterative method is called Modified New Iterative Method (MNIM) which is presented to enhance the validity domain of NIM techniques. The convergence analysis of the MNIM has also been illustrated. The simulation results show that the vaccination strategy can slow down the spread of the epidemic rapidly. Numerical results illustrate the excellent performance of the MNIM and show that the modified method is much more accurate than the NIM.

1. Introduction

The relation between the biology of infectious diseases, the propagation mechanism, and the mathematics used to explain them is a vital research undertaking. Infectious diseases frequently affect a significant number of populations across a wide geographical area. Experiments undertaken in laboratories are insufficient due to the vast differences in scale. Large-scale trials on infectious diseases in humans are not feasible or immoral [1]. Therefore, various mathematical models have been designed and utilized to obtain insight into the dynamics of transmission and control of the epidemic. To combat the risks of contagious epidemic diseases, it is critical to reducing the peak time of epidemic disease to slow down the inevitable dynamics and prepare for the subsequent epidemic wave [2].

In the past decades, the global scientific community has developed several mathematical models for the transmission and prevention of infectious diseases such as childhood disease [3], HIV/AIDS [4], dengue [5], tuberculosis [6], Ebola [7], COVID-19 [8,9,10,11,12,13] and many others. These models are very critical in diverse areas, such as policy making, emergency planning, risk management, prevention strategies, and the implementation of the enhancement of various health and economic aspects [14]. However, there is still a lack of numerical studies that can help the researcher to understand the disease transmission rate, infection rate, and recovery rate in the epidemic situation.

The susceptible-infected-recovered (SIR) model of epidemiological type, first proposed by Kermack and Mckendrick [15], is an epidemic model for the spread of infectious disease. In order to conduct the numerical study of the dynamics of such diseases, we considered the susceptible-infected-recovered (SIR) epidemiologic nonlinear mathematical model which is being created by O.D. Makinde [3] to track the temporal dynamics of a pediatric illness in the presence of a preventative vaccination strategy. In particular, compartmental relations are used to describe an infection state and a process of infection by dividing the population into compartments based on assumptions about the nature and rate of transfer from one compartment to another. The rate of leaving the compartment depends only on the current state of the epidemic [16].

Figure 1 demonstrates the SIR model of childhood infectious disease [17] contains a susceptible group, symbolized by $S$, an infected group $I$, and a recovered group $R$, which denotes a vaccinated, as well as recovered people with constant immunity. It assumes that the vaccine efficacy is 100% and the fraction of citizens vaccinated at birth, each year are denoted as $p$ (with $0<p<1$) and the remain populations $\left(1-p\right)$ are susceptible. A susceptible individual move into the infected group through communication with an infected individual by an average contact rate $b$. A person who has been infected, recovers at a rate of $c$ and enters recovered group. The recovered group also contain vaccinated individual. The differential equations for the SIR model follow the natural death rates $\mu$ and the birth rates $a$. Therefore, the size of population $N$ is not stable.

The rates of change of the three populations are determined by the following system of three coupled nonlinear ordinary differential equations [17] used in this analysis

$$S^{\prime }=\left(1-p\right)aN-\mu S\left(t\right)-b\frac{S\left(t\right)I\left(t\right)}{N},$$

(1)

$$I^{\prime }=b\frac{S\left(t\right)I\left(t\right)}{N}-\left(c+\mu \right)I\left(t\right),$$

(2)

$$R^{\prime }=paN+cI\left(t\right)-\mu R\left(t\right)$$

(3)

Here, all parameters $p, a, b, c, \mu$ are positive constants. The total population $N$ varies with respect to time $t$ (in year) as we assume that $a\ne \mu .$ To make sure that the overall population $N$ is always inconsistent, we add the equations contain in the system (1)–(3) to obtain

$$\frac{dN}{dt}=\left(a-\mu \right)N,$$

(4)

Each group’s size is proportional to the total population, rearranging the above system (1)–(3) by assuming $\frac{S}{N}=x, \frac{I}{N}=y, \frac{R}{N}=z$; we obtained the following expression:

$$x^{\prime }=\left(1-p\right)a-ax\left(t\right)-by\left(t\right)x\left(t\right),$$

(5)

$$y^{\prime }=by\left(t\right)x\left(t\right)-\left(c+a\right)y\left(t\right),$$

(6)

$$z^{\prime }=pa+cy\left(t\right)-az\left(t\right)$$

(7)

Therefore, this study 3 focused on solving the system (5)–(7) with initial conditions (Ics)—

$$x\left(0\right)=s_0,y\left(0\right)=i_0,z\left(0\right)=r_0$$

(8)

In our literature survey, we found that the approximate analytical solutions of the childhood SIR infectious disease model were presented first time by Makinde [3], using the Adomian decomposition method (ADM). Yildirim and Cherruault [18] derived the approximate analytical expression for the SIR model using the Homotopy perturbation method (HPM). Recently, Mungkasi [17] corrected the inaccurate solutions of Ghotbi et al. [19] by using successive approximation methods to the determined system. Srivastava and Günerhan [20] introduced the conformable fractional differential transform method (CFDTM) to calculate an approximate solution of the fractional-order susceptible-infected-recovered (SIR) epidemic model of childhood disease. The combined technique of the classical homotopy perturbation method and the Elzaki transform have been utilized by Jena et.al [21] to solve the childhood diseases model with fractional order.

In this research, we have revisited the model in order to provide a numerical framework via New Iterative Method (NIM) and Modified New Iterative Method (MNIM). The current study investigates the well-documented SIR infection model to test the validity domain of the proposed methods with straightforward computations. For this particular form of SIR model, the analytical solutions have been derived using a coalescence of Laplace transform and the New Iterative Method, which is called modified NIM. We also use the NIM scheme which is based on the iterative principle of characteristics, introduced by Daftardar-Gejji and Jafari [22]. Recently, the application of NIM can be found in the literature specifically dealing with Klein-Gordon equations [23], fractional Whitham-Broer-Kaup (WBK) system [24], Fisher’s equation [25], Cahn-Hillard equation [26], multidimensional wave equations [27], 1D Swift–Hohenberg equation [28], Falkner-Skan Equation [29], Jeffery-Hamel flow problem [30], Chemical kinetics equations [31,32], Lake pollution model [33] and many other problems. Stability analysis, Error bound, and Convergence analysis of NIM have been illustrated in Refs. [33,34].

2. Solution Procedure

2.1. New Iterative Method

In this section, a class of n-th-order differential equation has been considered in the form,

$$u^{\left(n\right)}\left(t\right)=f(t,u\left(t\right),u^{\prime }\left(t\right),u^{″}\left(t\right),\cdots u^{\left(n-1\right)}\left(t\right)),0\le t\le T,$$

(9)

Subject to the initial conditions

$$u\left(t_0\right)=u_0,u^{\prime }\left(t_0\right)=u_0^{\prime },u^{″}\left(t_0\right)=u_0^{″},\cdots u^{\left(n-1\right)}\left(t_0\right)=u_0{}^{n-1},$$

(10)

where $f$ is a time dependent, continuous linear/nonlinear function and $u_0,u_0^{\prime },u_0^{″},\cdots u_0{}^{n-1}$ are prescribed. Now the operator form of the above equation can be written as,

$$Lu=R\left(u\right)+N\left(u\right)+g\left(t\right),$$

(11)

where $f\left(t,u\left(t\right),u^{\prime }\left(t\right),u^{″}\left(t\right),\cdots u^{\left(n-1\right)}\left(t\right)\right)=Ru+Nu+g\left(t\right),$ and $L\equiv \frac{d^n}{d{t}^n}$ is chosen to be highest order derivative, which is easily invertible, $Ru$ the remaining linear differential operator whose order is less than $L$; $Nu$ represents the nonlinear terms and $g$ is the given inhomogeneous source term. The Equation (11) can now be rewritten as:

$$Lu=g+Ru+Nu,$$

(12)

Applying the inverse operator $L^{-1}$ to both sides of (12),

$$u=f\left(t\right)+L^{-1}Ru+L^{-1}Nu,$$

(13)

where $f\left(t\right)$ is the term arising from integrating the source term $g\left(t\right)$.

Let, the solution of Equation (13) for $u$ may be written in the series:

$$u\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{u}_i=u_0+u_1+u_2+\dots$$

(14)

The nonlinear term $N$ can be find out as

$$H_0=N\left(u_0\right) \mathrm{and} H_m=N\left({\displaystyle \sum }_{j=0}^m{u}_j\right)-N\left({\displaystyle \sum }_{j=0}^{m-1}{u}_j\right).$$

(15)

Thus, according to NIM,

$$N\left({{\displaystyle \sum }}_{i=0}^{\infty }{u}_i\right)=H_0+H_1+H_2+\dots =N\left(u_0\right)+\left(N\left(u_0+u_1\right)-N\left(u_0\right)\right)+\left(N\left(u_0+u_1+u_2\right)-N\left(u_0+u_1\right)\right)+\dots$$

(16)

Now, we define the recurrence relation as follows:

$$u_0=L^{-1}g\left(t\right)+\phi =f\left(t\right),$$

(17)

$$u_1=L^{-1}R\left(u_0\right)+L^{-1}\left(H_0\right),$$

(18)

$$u_{m+1}=L^{-1}R\left(u_m\right)+L^{-1}\left(H_m\right), m\ge 1$$

(19)

Since R is linear then,

$$\left({\displaystyle \sum }_{i=0}^{m}R{u}_i\right)=R\left({\displaystyle \sum }_{i=0}^{\infty }{u}_i\right)$$

(20)

Then series solution becomes

$${\displaystyle \sum }_{i=1}^{\infty }{u}_i=f\left(t\right)+L^{-1}R\left({\displaystyle \sum }_{i=0}^{\infty }{u}_i\right)+L^{-1}N\left({\displaystyle \sum }_{i=0}^{\infty }{u}_i\right)$$

(21)

Now, the n-term solution can be written as series form as:

$$u={\displaystyle \sum }_{i=0}^{n-1}{u}_i=u_0+u_1+u_2+\dots +u_{n-1}.$$

(22)

2.2. Modified New Iterative Method

We have proposed an improvement of NIM based on the composition of Laplace transform with the New Iterative Method to solve system of differential equations.

Let us consider a system of differential equations in the operator form:

$$D\left({\mu }_1\right)+R_1\left({\mu }_1\right)+N_1\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)=g_1\left(t\right),$$

(23)

$$D\left({\mu }_2\right)+R_2\left({\mu }_2\right)+N_2\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)=g_2\left(t\right),$$

(24)

$$\begin{array}{c}⋮\\ D\left({\mu }_n\right)+R_n\left({\mu }_n\right)+N_n\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)=g_n\left(t\right),\end{array}$$

(25)

Subjected to the initial conditions.

$${\mu }_1\left(0\right)=a_1,{\mu }_2\left(0\right)=a_2,...,{\mu }_n\left(0\right)=a_n,$$

(26)

where $D$ is a invertible linear differential operator, and $R_1,R_2,...,R_n$ are remaining linear operators order less than $D$ and $N_1,N_2,...,N_n$ are nonlinear operators and $g_1,g_2,...,g_n$ are inhomogeneous terms. For the first order $D\equiv \frac{d}{dt}$, second order $D\equiv \frac{d^2}{d{t}^2}$ and so on.

The technique consists first of applying Laplace transformation (which is denoted by $\Im$) to both sides of systems (23)–(25), hence

$$\Im \left\{D\left({\mu }_1\right)\right\}+\Im \left\{R_1\left({\mu }_1\right)\right\}+\Im \left\{N_1\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\}=\Im \left\{g_1\left(t\right)\right\},$$

(27)

$$\Im \left\{D\left({\mu }_2\right)\right\}+\Im \left\{R_2\left({\mu }_2\right)\right\}+\Im \left\{N_2\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\}=\Im \left\{g_2\left(t\right)\right\},$$

(28)

$$\begin{array}{c}⋮\\ \Im \left\{D\left({\mu }_n\right)\right\}+\Im \left\{R_n\left({\mu }_n\right)\right\}+\Im \left\{N_n\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\}=\Im \left\{g_n\left(t\right)\right\}\end{array}$$

(29)

Applying the formulas for Laplace transforms, we obtain

$$s\Im \left({\mu }_1\right)-{\mu }_1\left(0\right)=\Im \left\{g_1\left(t\right)\right\}-\Im \left\{R_1\left({\mu }_1\right)\right\}-\Im \left\{N_1\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\},$$

(30)

$$s\Im \left({\mu }_2\right)-{\mu }_2\left(0\right)=\Im \left\{g_2\left(t\right)\right\}-\Im \left\{R_2\left({\mu }_2\right)\right\}-\Im \left\{N_2\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\},$$

(31)

$$\begin{array}{c}⋮\\ s\Im \left({\mu }_n\right)-{\mu }_n\left(0\right)=\Im \left\{g_n\left(t\right)\right\}-\Im \left\{R_n\left({\mu }_n\right)\right\}-\Im \left\{N_n\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\},\end{array}$$

(32)

where ‘$s$’ is called a Laplace domain function.

Using the initial conditions (26), we have,

$$\Im \left({\mu }_1\right)=\frac{a_1}{s}+\frac{\Im \left\{g_1\left(t\right)\right\}}{s}-\frac{1}{s}\Im \left\{R_1\left({\mu }_1\right)\right\}-\frac{1}{s}\Im \left\{N_1\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\},$$

(33)

$$\Im \left({\mu }_2\right)=\frac{a_2}{s}+\frac{\Im \left\{g_2\left(t\right)\right\}}{s}-\frac{1}{s}\Im \left\{R_2\left({\mu }_2\right)\right\}-\frac{1}{s}\Im \left\{N_2\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\},$$

(34)

$$\begin{array}{c}⋮\\ \Im \left({\mu }_n\right)=\frac{a_n}{s}+\frac{\Im \left\{g_n\left(t\right)\right\}}{s}-\frac{1}{s}\Im \left\{R_n\left({\mu }_n\right)\right\}-\frac{1}{s}\Im \left\{N_n\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\},\end{array}$$

(35)

Applying the inverse Laplace transform to the equations in (33)–(35), we get

$${\mu }_1={\Im }^{-1}\left[\frac{a_1}{s}+\frac{\Im \left\{g_1\left(t\right)\right\}}{s}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left\{R_1\left({\mu }_1\right)\right\}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left\{N_1\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\}\right],$$

(36)

$${\mu }_2={\Im }^{-1}\left[\frac{a_2}{s}+\frac{\Im \left\{g_2\left(t\right)\right\}}{s}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left\{R_2\left({\mu }_2\right)\right\}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left\{N_2\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\}\right],$$

(37)

$$\begin{array}{c}⋮\\ {\mu }_n={\Im }^{-1}\left[\frac{a_2}{s}+\frac{\Im \left\{g_n\left(t\right)\right\}}{s}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left\{R_n\left({\mu }_n\right)\right\}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left\{N_n\left({\mu }_1,{\mu }_2,...,{\mu }_n\right)\right\}\right]\end{array}$$

(38)

Let, the approximate the solutions ${\mu }_1,{\mu }_2,...,{\mu }_n$ of system (36)–(38) can be expressed as

$${\mu }_1={\mu }_{1,0}+{\mu }_{1,1}+\cdots ={\displaystyle \sum }_{i=0}^{\infty }{\mu }_{1,i},$$

(39)

$${\mu }_2={\mu }_{2,0}+{\mu }_{2,1}+\cdots ={\displaystyle \sum }_{i=0}^{\infty }{\mu }_{2,i},$$

(40)

$$\begin{array}{c}⋮\\ {\mu }_n={\mu }_{n,0}+{\mu }_{n,1}+\cdots ={\displaystyle \sum }_{i=0}^{\infty }{\mu }_{n,i}.\end{array}$$

(41)

For the NIM, the nonlinear operators can be decomposed by,

$$G_0=N_1\left({\mu }_{1,0}\right) \mathrm{and} G_m=N_1\left({{\displaystyle \sum }}_{i=0}^m{\mu }_{1,i}\right)-N_1\left({{\displaystyle \sum }}_{i=0}^{m-1}{\mu }_{1,i}\right),m\ge 1,$$

(42)

$$H_0=N_2\left({\mu }_{2,0}\right) \mathrm{and} H_m=N_2\left({{\displaystyle \sum }}_{i=0}^m{\mu }_{2,i}\right)-N_2\left({{\displaystyle \sum }}_{i=0}^{m-1}{\mu }_{2,i}\right),m\ge 1,$$

(43)

$$\begin{array}{c}⋮\\ I_0=N_n\left({\mu }_{n,0}\right) \mathrm{and} I_m=N_n\left({{\displaystyle \sum }}_{i=0}^m{\mu }_{n,i}\right)-N_n\left({{\displaystyle \sum }}_{i=0}^{m-1}{\mu }_{n,i}\right),m\ge 1\end{array}$$

(44)

Thus, the nonlinear operators can be expressed as:

$$N_1\left({{\displaystyle \sum }}_{i=0}^{\infty }{\mu }_{1,i}\right)=N_1\left({\mu }_{1,0}\right)+\left\{N_1\left({\mu }_{1,0}+{\mu }_{1,1}\right)-N_1\left({\mu }_{1,0}\right)\right\}+\left\{N_1\left({\mu }_{1,0}+{\mu }_{1,1}+{\mu }_{1,2}\right)-N_1\left({\mu }_{1,0}+{\mu }_{1,1}\right)\right\}+\cdots ,$$

(45)

$$N_2\left({{\displaystyle \sum }}_{i=0}^{\infty }{\mu }_{2,i}\right)=N_2\left({\mu }_{2,0}\right)+\left\{N_2\left({\mu }_{2,0}+{\mu }_{2,1}\right)-N_2\left({\mu }_{2,0}\right)\right\}+\left\{N_2\left({\mu }_{2,0}+{\mu }_{2,1}+{\mu }_{2,2}\right)-N_2\left({\mu }_{2,0}+{\mu }_{2,1}\right)\right\}+\cdots ,$$

(46)

$$\begin{array}{c}⋮\\ N_n\left({{\displaystyle \sum }}_{i=0}^{\infty }{\mu }_{n,i}\right)=N_n\left({\mu }_{n,0}\right)+\left\{N_n\left({\mu }_{n,0}+{\mu }_{n,1}\right)-N_n\left({\mu }_{n,0}\right)\right\}+\left\{N_n\left({\mu }_{n,0}+{\mu }_{n,1}+{\mu }_{n,2}\right)-N_n\left({\mu }_{n,0}+{\mu }_{n,1}\right)\right\}+\cdots \end{array}$$

(47)

Since $R_1,R_2,...,R_n$ are linear.

$${\displaystyle \sum }_{i=0}^{\infty }{R}_1\left({\mu }_{1,i}\right)=R_1\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{1,i}\right),$$

(48)

$${\displaystyle \sum }_{i=0}^{\infty }{R}_2\left({\mu }_{2,i}\right)=R_2\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{2,i}\right),$$

(49)

$$\begin{array}{c}⋮\\ {\displaystyle \sum _{i=0}^{\infty }}{R}_n\left({\mu }_{n,i}\right)=R_n\left({\displaystyle \sum _{i=0}^{\infty }}{\mu }_{n,i}\right)\end{array}$$

(50)

The NIM admits the use of the recursive relations by following way:

$${\mu }_{1,0}={\Im }^{-1}\left[\frac{\Im \left\{g_1\left(t\right)\right\}}{s}\right]+a_1=f_1\left(t\right)+a_1,$$

(51)

where $f_1\left(t\right)$ is the term arising after inverse Laplace transformation of the source term $\frac{\Im \left\{g_1\left(t\right)\right\}}{s},$ all of which are assumed to be prescribed.

$${\mu }_{1,1}=-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_1\left({\mu }_{1,0}\right)\right\}\right]-\frac{1}{s}{\Im }^{-1}\left[\Im \left(G_0\right)\right],$$

(52)

$$\begin{array}{c}⋮\\ {\mu }_{1,m+1}=-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_1\left({\mu }_{1,m}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left(G_m\right)\right],m=1,2,....\end{array}$$

(53)

Thus,

$$\therefore {\mu }_1={\displaystyle \sum }_{i=0}^{\infty }{\mu }_{1,i}=f_1+{\displaystyle \sum }_{i=1}^{\infty }{\mu }_{1,i}=f_1-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_1\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{1,i}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{N_1\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{1,i}\right)\right\}\right]$$

(54)

By the similar manner, we have

$${\mu }_{2,0}={\Im }^{-1}\left[\frac{\Im \left\{g_2\left(t\right)\right\}}{s}\right]+a_2=f_2\left(t\right)+a_2,$$

(55)

$${\mu }_{2,1}=-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_2\left({\mu }_{2,0}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left(H_0\right)\right],$$

(56)

$$\begin{array}{c}⋮\\ {\mu }_{2,m+1}=-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_2\left({\mu }_{2,m}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left(H_m\right)\right],m=1,2,...,\end{array}$$

(57)

$$\therefore {\mu }_2={\displaystyle \sum }_{i=0}^{\infty }{\mu }_{2,i}=f_2+{\displaystyle \sum }_{i=1}^{\infty }{\mu }_{2,i}=f_2-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_2\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{2,i}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{N_2\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{2,i}\right)\right\}\right],$$

(58)

and

$${\mu }_{n,0}={\Im }^{-1}\left[\frac{\Im \left\{g_n\left(t\right)\right\}}{s}\right]+a_n=f_n\left(t\right)+a_n,$$

(59)

$${\mu }_{n,1}=-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_n\left({\mu }_{n,0}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{I_0\right\}\right],$$

(60)

$$\begin{array}{c}⋮\\ {\mu }_{n,m+1}=-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_n\left({\mu }_{n,m}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{I_m\right\}\right],m=1,2,...,\end{array}$$

(61)

$$\therefore {\mu }_n={\displaystyle \sum }_{i=0}^{\infty }{\mu }_{n,i}=f_n+{\displaystyle \sum }_{i=1}^{\infty }{\mu }_{n,i}=f_n-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{R_n\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{n,i}\right)\right\}\right]-{\Im }^{-1}\left[\frac{1}{s}\Im \left\{N_n\left({\displaystyle \sum }_{i=0}^{\infty }{\mu }_{n,i}\right)\right\}\right]$$

(62)

The $k$-terms approximate solutions are given by the following form:

$${\mu }_1={\displaystyle \sum }_{i=0}^{k-1}{\mu }_{1,i}={\mu }_{1,0}+{\mu }_{1,1}+{\mu }_{1,2}+\cdots +{\mu }_{1,k-1},$$

(63)

$${\mu }_2={\displaystyle \sum }_{i=0}^{k-1}{\mu }_{2,i}={\mu }_{2,0}+{\mu }_{2,1}+{\mu }_{2,2}+\cdots +{\mu }_{2,k-1},$$

(64)

$$\begin{array}{c}⋮\\ {\mu }_n={\displaystyle \sum _{i=0}^{k-1}}{\mu }_{n,i}={\mu }_{n,0}+{\mu }_{n,1}+{\mu }_{n,2}+\cdots +{\mu }_{n,k-1}.\end{array}$$

(65)

2.3. Convergence Analysis of MNIM

Let, ${\rm X},\Upsilon$ and $\kappa$ be the elements in a Banach space $B$ and $N$ is nonlinear contraction from $B\to B$ such that $\Vert \kappa \Vert =\Vert N\left({\rm X}\right)-N\left(\Upsilon \right)\Vert \le \sigma \Vert {\rm X}-\Upsilon \Vert , 0<\sigma <1$. Then we can prove that $\Vert {\kappa }_{r+1}\Vert \le {\sigma }^{r+1}\Vert {\kappa }_0\Vert ,r=0,1,...$ by the principle of Banach fixed point theorem.

Let, ${\kappa }_0$ is the initial approximation and $\rho$ is a known function in Banach Space $B$.

$${\kappa }_0\left(t\right)=\rho \left(t\right),$$

(66)

By the decomposition principle of nonlinear terms in MNIM, we can write

$$\Vert {\kappa }_1\Vert =\Vert N\left({\kappa }_0\right)\Vert \le \sigma \Vert {\kappa }_0\Vert ,$$

(67)

$$\Vert {\kappa }_2\Vert =\Vert N\left({\kappa }_0+{\kappa }_1\right)-N\left({\kappa }_0\right)\Vert \le \sigma \Vert {\kappa }_1\Vert \le {\sigma }^2\Vert {\kappa }_0\Vert ,$$

(68)

$$\Vert {\kappa }_3\Vert =\Vert N\left({\kappa }_0+{\kappa }_1+{\kappa }_2\right)-N\left({\kappa }_0+{\kappa }_1\right)\Vert \le \sigma \Vert {\kappa }_2\Vert \le {\sigma }^2\Vert {\kappa }_1\Vert \le {\sigma }^3\Vert {\kappa }_0\Vert ,$$

(69)

$$\begin{array}{c}⋮\\ \Vert {\kappa }_{r+1}\Vert =\Vert N({\displaystyle \sum }_{j=0}^r{\kappa }_j)-N({\displaystyle \sum }_{j=0}^{r-1}{\kappa }_j)\Vert \le \sigma \Vert {\kappa }_r\Vert \le {\sigma }^{r+1}\Vert {\kappa }_0\Vert ,r=0,1,2,....\end{array}$$

(70)

Hence, the solutions of MNIM procedure converges via Banach fixed point theorem.

3. Application

3.1. Solutions by NIM

To check the validity of NIM solutions we have solved the problem by NIM method. Let us take the integration on the system (5)–(7), we obtain:

$$x\left(t\right)=s_0+\left(1-p\right)at-\left({{\displaystyle \int }}_0^{t}ax\left(t\right)dt\right)-\left({{\displaystyle \int }}_0^{t}by\left(t\right)x\left(t\right)dt\right),$$

(71)

$$y\left(t\right)=i_0+\left({{\displaystyle \int }}_0^{t}by\left(t\right)x\left(t\right)dt\right)-\left({{\displaystyle \int }}_0^t\left(c+a\right)y\left(t\right)dt\right),$$

(72)

$$z\left(t\right)=r_0+pat+\left({{\displaystyle \int }}_0^{t}cy\left(t\right)dt\right)-\left({{\displaystyle \int }}_0^{t}az\left(t\right)dt\right)$$

(73)

In view of Equations (17)–(19) and according to the solution procedure of NIM for system of differential equations, we find the approximations as follows:

$$x_0\left(t\right)=s_0+\left(1-p\right)at,$$

(74)

$$y_0\left(t\right)=i_0,$$

(75)

$$z_0\left(t\right)=pat+r_0,$$

(76)

$$x_1\left(t\right)=-\frac{\left(1-p\right)a^2{t}^2}{2}-a{s}_{0}t-\frac{b{i}_0\left(1-p\right)a{t}^2}{2}-b{i}_0{s}_{0}t,$$

(77)

$$y_1\left(t\right)=\frac{b{i}_0\left(1-p\right)a{t}^2}{2}+b{i}_0{s}_{0}t-\left(c+a\right)i_{0}t,$$

(78)

$$z_1\left(t\right)=c{i}_{0}t-\frac{1}{2}p{a}^2{t}^2-a{r}_{0}t,$$

(79)

Similarly, we can obtain the other approximations as well. Therefore, the approximate precisions are:

$$x\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{x}_i\left(t\right),$$

(80)

$$y\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{y}_i\left(t\right),$$

(81)

$$z\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{z}_i\left(t\right).$$

(82)

3.2. Solutions by MNIM

By applying Laplace Transform on the system of Equations (5)–(7), and using ICs (8), we obtain:

$$\Im \left(x\left(t\right),t,\epsilon \right)=\frac{1}{\epsilon +a}\left(-\frac{a\left(p-1\right)}{\epsilon }+s_0-b\Im \left(y\left(t\right) x\left(t\right),t,\epsilon \right)\right),$$

(83)

$$\Im \left(y\left(t\right),t,\epsilon \right)=\frac{1}{\epsilon +a+c}\left(i_0+b\Im \left(y\left(t\right)x\left(t\right),t,\epsilon \right)\right),$$

(84)

$$\Im \left(z\left(t\right),t,\epsilon \right)=\frac{1}{\epsilon +a}\left(\frac{pa}{\epsilon }+r_0+c\Im \left(y\left(t\right) ,t,\epsilon \right)\right),$$

(85)

Taking inverse Laplace transform operator on (83)–(85), we get:

$$x\left(t\right)={\Im }^{-1}\left\{\frac{1}{\epsilon +a}\left(-\frac{a\left(p-1\right)}{\epsilon }+s_0-b\Im \left(y\left(t\right)x\left(t\right),t,\epsilon \right)\right)\right\},$$

(86)

$$y\left(t\right)={\Im }^{-1}\left\{\frac{1}{\epsilon +a+c}\left(i_0+b\Im \left(y\left(t\right)x\left(t\right),t,\epsilon \right)\right)\right\},$$

(87)

$$z\left(t\right)={\Im }^{-1}\left\{\frac{1}{\epsilon +a}\left(\frac{pa}{\epsilon }+r_0+c\Im \left(y\left(t\right) ,t,\epsilon \right)\right)\right\},$$

(88)

In view of (51)–(62), we can obtain the following iterations by using MNIM,

$$x_0\left(t\right)={\Im }^{-1}\left\{\frac{1}{\epsilon +a}\left(-\frac{a\left(p-1\right)}{\epsilon }+s_0\right)\right\}=1-p+\left(p+s_0-1\right)e^{-at},$$

(89)

$$y_0\left(t\right)={\Im }^{-1}\left\{\frac{i_0}{\epsilon +a+c}\right\}=i_0{e}^{-\left(c+a\right)t},$$

(90)

$$z_0\left(t\right)={\Im }^{-1}\left\{\frac{\frac{pa}{\epsilon }+r_0}{\epsilon +a}\right\}=p+\left(-p+r_0\right)e^{-at},$$

(91)

$$x_1\left(t\right)=\frac{\left(e^{-\left(c+2a\right)t}c\left(p+s_0-1\right)-\left(p-1\right)e^{-\left(c+a\right)t}\left(c+a\right)+e^{-at}\left(pa-c{s}_0-a\right)\right)i_{0}b}{c\left(c+a\right)},$$

(92)

$$y_1\left(t\right)=\frac{\left(e^{-\left(c+2a\right)t}\left(-p-s_0+1\right)+\left(-apt+at+p+s_0-1\right)e^{-\left(c+a\right)t}\right)i_{0}b}{a},$$

(93)

$$z_1\left(t\right)=i_0\left(-e^{-\left(c+a\right)t}+e^{-at}\right),$$

(94)

and so on. In the similar approach the further components can be obtained. Therefore, the approximate precisions are

$$x\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{x}_i\left(t\right),$$

(95)

$$y\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{y}_i\left(t\right),$$

(96)

$$z\left(t\right)={\displaystyle \sum }_{i=0}^{\infty }{z}_i\left(t\right).$$

(97)

4. Results and Discussions

We take the initial conditions $s_0=0.8,i_0=0.2,r_0=0$ and parameters values $p=0.9,a=0.4,b=0.8,$ $c=0.03$; a similar scenario of ref. [3]. All our calculations as well as our graphs are carried out by Maple 2020. We set the continuous constant step size $h=0.001$ for the RK4 in the maple software.

4.1. NIM Solutions

Therefore, 5-iterations of NIM solutions are:

$$\begin{array}{c}x\left(t\right)={\displaystyle \sum _{i=0}^5}{x}_i\left(t\right)\\ =6.736115226612228\times {10}^{-52}{t}^{47}+3.922133162637945\times {10}^{-49}{t}^{46}\\ +1.073655848446378\times {10}^{-46}{t}^{45}+1.844655274257630\times {10}^{-44}{t}^{44}\\ +2.242032611955568\times {10}^{-42}{t}^{43}+2.060673733084379\times {10}^{-40}{t}^{42}\\ +1.494716735333350\times {10}^{-38}{t}^{41}+8.809612654525537\times {10}^{-37}{t}^{40}\\ +4.305357070507673\times {10}^{-35}{t}^{39}\\ +1.768963027457917\times {10}^{-33}{t}^{38}+6.164360273888642\times {10}^{-32}{t}^{37}\\ +1.829859631340267\times {10}^{-30}{t}^{36}+4.627443481597180\times {10}^{-29}{t}^{35}\\ +9.920911834830423\times {10}^{-28}{t}^{34}+1.782386359157935\times {10}^{-26}{t}^{33}\\ +2.622401070812103\times {10}^{-25}{t}^{32}+3.013741998638641\times {10}^{-24}{t}^{31}\\ +2.397281814414504\times {10}^{-23}{t}^{30}+7.028239936204403\times {10}^{-23}{t}^{29}\\ -1.230670827410604\times {10}^{-21}{t}^{28}-2.070457639465569\times {10}^{-20}{t}^{27}\\ -1.239142909703323\times {10}^{-19}{t}^{26}+2.949835075276537\times {10}^{-19}{t}^{25}\\ +3.725341255626546\times {10}^{-17}{t}^{23}-3.435588549817450\times {10}^{-16}{t}^{22}\\ +9.557233694959134\times {10}^{-18}{t}^{24}-2.904383929261819\times {10}^{-15}{t}^{21}\\ +1.169109476870975\times {10}^{-14}{t}^{20}+1.575620033256645\times {10}^{-13}{t}^{19}\\ -4.676372541642090\times {10}^{-13}{t}^{18}-6.907013200303840\times {10}^{-12}{t}^{17}\\ +2.695714904867200\times {10}^{-11}{t}^{16}+1.988585322443935\times {10}^{-10}{t}^{15}\\ -1.490922789811607\times {10}^{-9}{t}^{14}-6.672613505299390\times {10}^{-10}{t}^{13}\\ +3.610221503703013\times {10}^{-8}{t}^{12}-5.718040776900796\times {10}^{-8}{t}^{11}\\ -7.848058588316446\times {10}^{-7}{t}^{10}+4.260175583211709\times {10}^{-6}{t}^9\\ -1.110385716627811\times {10}^{-6}{t}^8-0.00005754643154102856t^7\\ +0.0001517085007466667t^6+0.0002838500160000001t^5\\ -0.001811776t^4-0.008224t^3+0.1008t^2-0.408t+0.8,\end{array}$$

(98)

$$\begin{array}{c}y\left(t\right)={\displaystyle \sum _{i=0}^5}{y}_i\left(t\right)\\ =6.736115226612228\times {10}^{-52}{t}^{47}-3.922133162637945\times {10}^{-49}{t}^{46}\\ -1.073655848446378\times {10}^{-46}{t}^{45}-1.844655274257630\times {10}^{-44}{t}^{44}\\ -2.242032611955568\times {10}^{-42}{t}^{43}-2.060673733084379\times {10}^{-40}{t}^{42}\\ -1.494716735333350\times {10}^{-38}{t}^{41}-8.809612654525537\times {10}^{-37}{t}^{40}\\ -4.305357070507673\times {10}^{-35}{t}^{39}\\ -1.768963027457917\times {10}^{-33}{t}^{38}-6.164360273888642\times {10}^{-32}{t}^{37}\\ -1.829859631340267\times {10}^{-30}{t}^{36}-4.627443481597180\times {10}^{-29}{t}^{35}\\ -9.920911834830423\times {10}^{-28}{t}^{34}-1.782386359157935\times {10}^{-26}{t}^{33}\\ -2.622401070812103\times {10}^{-25}{t}^{32}-3.013741998638641\times {10}^{-24}{t}^{31}\\ -2.397281814414504\times {10}^{-23}{t}^{30}-7.028239936204403\times {10}^{-23}{t}^{29}\\ +1.230670827410604\times {10}^{-21}{t}^{28}+2.070457639465569\times {10}^{-20}{t}^{27}\\ +1.239142909703323\times {10}^{-19}{t}^{26}-2.949835075276537\times {10}^{-19}{t}^{25}\\ -9.557233694710467\times {10}^{-18}{t}^{24}-3.725341248233134\times {10}^{-17}{t}^{23}\\ +3.435588646687410\times {10}^{-16}{t}^{22}+2.904384680530290\times {10}^{-15}{t}^{21}\\ -1.169105545529128\times {10}^{-14}{t}^{20}-1.575605121898856\times {10}^{-13}{t}^{19}\\ +4.676793838257692\times {10}^{-13}{t}^{18}+6.907896025018594\times {10}^{-12}{t}^{17}\\ -2.694385690830112\times {10}^{-11}{t}^{16}-1.987319350530539\times {10}^{-10}{t}^{15}\\ +1.491344335393742\times {10}^{-9}{t}^{14}+6.609946732092200\times {10}^{-10}{t}^{13}\\ -3.617958560529853\times {10}^{-8}{t}^{12}+5.712957034873109\times {10}^{-8}{t}^{11}\\ +7.885699621920298\times {10}^{-7}{t}^{10}-4.260056651414198\times {10}^{-6}{t}^9\\ +9.012201347078117\times {10}^{-7}{t}^8+0.00005827143817874285t^7\\ -0.0001487855181688889t^6-0.0003762609485t^5\\ +0.002759918749999999t^4-0.0011697t^3-0.02823t^2\\ +0.042t+0.2;\end{array}$$

(99)

$$\begin{array}{c}z\left(t\right)={\displaystyle \sum _{i=0}^5}{z}_i\left(t\right)\\ =-2.486673133813368\times {10}^{-28}{t}^{24}-7.393412085355813\times {10}^{-26}{t}^{23}\\ -9.686995915539818\times {10}^{-24}{t}^{22}-7.512684713956981\times {10}^{-22}{t}^{21}\\ -3.931341846821115\times {10}^{-20}{t}^{20}-1.491135778969304\times {10}^{-18}{t}^{19}\\ -4.212966156027532\times {10}^{-17}{t}^{18}-8.828247147539671\times {10}^{-16}{t}^{17}\\ -1.329214037088414\times {10}^{-14}{t}^{16}-1.265971913395598\times {10}^{-13}{t}^{15}\\ -4.215455821345592\times {10}^{-13}{t}^{14}+6.266677320718978\times {10}^{-12}{t}^{13}\\ +7.737056826840045\times {10}^{-11}{t}^{12}+5.083742027686752\times {10}^{-11}{t}^{11}\\ -3.764103360385219\times {10}^{-9}{t}^{10}-1.18931797511109\times {10}^{-10}{t}^9\\ +2.091655819199999\times {10}^{-7}{t}^8-7.250066377142856\times {10}^{-7}{t}^7\\ -8.611871466666666\times {10}^{-6}{t}^6+0.00009241093249999999t^5\\ -0.00094814275t^4+0.0093937t^3-0.07257t^2+0.366t\end{array}$$

(100)

4.2. MNIM Solutions

The 3-iterations MNIM solutions are:

$$\begin{array}{c}6x\left(t\right)={\displaystyle \sum _{i=0}^3}{x}_i\left(t\right)\\ =0.1-0.0006017646192844911e^{-3.32t}+0.03925876631627679e^{-1.63t}\\ +0.0007765131412058982e^{-2.89t}-0.0002581915613754199e^{-2.46t}\\ -8.154807242157790e^{-0.4t}+1.535711286278180\times {10}^{-23}{e}^{-1.66t}\\ \times (3.233242835823126\times {10}^{17}{t}^2+3.406539138328771\times {10}^{20}t\\ -4.593848682004663\times {10}^{21})-5.031003559435018\times {10}^{-22}{e}^{-0.83 t}\\ \times (5.439324414583916\times {10}^{18}{t}^2+3.505486748642923\times {10}^{20}t\\ +5.176936742877226\times {10}^{21})+4.232804232804233\times {10}^{-20}{e}^{-2.92t}\\ \times \left(2.774146215816360\times {10}^{15}t-1.179728723726773\times {10}^{16}\right)\\ -8.597701458250771\times {10}^{-25}{e}^{-2.52t}\times (6.834618772122878\times {10}^{18}{t}^2\\ -6.887323759231177\times {10}^{20}t-2.054789750054483\times {10}^{22})\\ -9.360006622204685\times {10}^{-26}{e}^{-1.26t}\times (1.775957335388708\times {10}^{19}{t}^3\\ +7.601863095513286\times {10}^{20}{t}^2-4.338292456048305\times {10}^{23}t\\ -1.563308625518170\times {10}^{25})-2.012401423774007\times {10}^{-22}{e}^{-2.12t}\\ \times (2.800264785138102\times {10}^{17}{t}^2+8.054649373251765\times {10}^{18}t\\ -3.731099912706602\times {10}^{19})-3.205609817180065\times {10}^{-22}{e}^{-1.72t}\\ \times (4.454474776080082\times {10}^{17}{t}^2+2.940605751880004\times {10}^{19}t\\ +4.855607088780105\times {10}^{20})+7.623376607907109\times {10}^{-23}{e}^{-0.86t}\\ \times (6.229503999994038\times {10}^{17}{t}^3+8.089326933329536\times {10}^{19}{t}^2\\ +3.188740492482025\times {10}^{21}t+3.781029478577795\times {10}^{22})\\ +1.481481481481481\times {10}^{-19}{e}^{-0.43t}\times (1.151999999999448\times {10}^{16}{t}^2\\ +1.459199999999788\times {10}^{18}t+5.799796492739982\times {10}^{19})\\ -1.831459902474760\times {10}^{-22}{e}^{-2.49t}\times (9.925997520612579\times {10}^{17}t\\ +2.746991500144911\times {10}^{19})-8.287048441320135\times {10}^{-26}{e}^{-2.09t}\\ \times (3.874906263024820\times {10}^{19}{t}^2-2.895379542793813\times {10}^{21}t\\ -4.195395911283663\times {10}^{23})+1.444445466389612\times {10}^{-27}{e}^{-1.69t}\\ \times (1.065222144015543\times {10}^{21}{t}^3+2.599021250303810\times {10}^{22}{t}^2\\ -5.333298360890297\times {10}^{23}t+7.995758255225330\times {10}^{25})\\ +4.250193528577767\times {10}^{-27}{e}^{-1.29t}\times (2.008732537975042\times {10}^{21}{t}^3\\ +4.866530225436270\times {10}^{22}{t}^2-8.085119609709582\times {10}^{24}t\\ -2.498760580734028\times {10}^{26})-1.161271592393671\times {10}^{-20}{e}^{-1.23t}\\ \times \left(5.725610654715977\times {10}^{17}t+3.019834074926216\times {10}^{19}\right)\\ +9.290172739149369\times {10}^{-22}{e}^{-2.06t}\times (7.745644989776837\times {10}^{16}t\\ -5.254962364229557\times {10}^{19}),\end{array}$$

(101)

$$\begin{array}{c}y\left(t\right)={\displaystyle \sum _{i=0}^3}{y}_i\left(t\right)\\ =0.0006080113108341571e^{-3.32t}+0.0002620072002134803e^{-2.46t}\\ -0.0007859828136596286e^{-2.89t}-0.04024023547418371e^{-1.63t}\\ +0.2e^{-0.43 t}+2.146561518692274\times {10}^{-23}{e}^{-1.72t}\\ \times (6.806879223943861\times {10}^{18}{t}^2+4.495935696722963\times {10}^{20}t\\ +7.427951840196584\times {10}^{21})-1.248000882960625\times {10}^{-23}{e}^{-0.86t}\\ \times (4.070758399996104\times {10}^{18}{t}^3+5.304608426664167\times {10}^{20}{t}^2\\ +2.100624173108962\times {10}^{22}t+2.506300632140585\times {10}^{23})\\ +5.453252023119212\times {10}^{-73}{e}^{-0.43t}(3.129631015457178\times {10}^{67}{t}^3\\ +2.816667913913444\times {10}^{69}{t}^2+7.626816999579977\times {10}^{70}t\\ +5.099515294485344\times {10}^{71})+3.770383636535017\times {10}^{-22}{e}^{-2.49t}\\ \times \left(4.891759543651175\times {10}^{17}t+1.354121357567373\times {10}^{19}\right)\\ +3.497806001185756\times {10}^{-25}{e}^{-2.09t}\times (9.346391669051500\times {10}^{18}{t}^2\\ -6.981744507856679\times {10}^{20}t-1.012015721679164\times {10}^{23})\\ -1.269602584555373\times {10}^{-26}{e}^{-1.69t}(1.240774054577326\times {10}^{20}{t}^3\\ +3.034218295207262\times {10}^{21}{t}^2-6.199976787711437\times {10}^{22}t\\ +9.312434748632805\times {10}^{24})-3.900002759251952\times {10}^{-26}{e}^{-1.29t}\\ \times (2.265465543829327\times {10}^{20}{t}^3+5.515152550771495\times {10}^{21}{t}^2\\ -9.113544300624442\times {10}^{23}t-2.821640082696647\times {10}^{25})\\ +2.5\times {10}^{-19}{e}^{-1.23t}\times (2.759330436007700\times {10}^{16}t\\ +1.456588408293963\times {10}^{18})-4.817644623433325\times {10}^{-22}{e}^{-2.06t}\\ \times \left(1.521132690763403\times {10}^{17}t-1.031981767729906\times {10}^{20}\right)\\ -4.127108492976661\times {10}^{-24}{e}^{-1.66t}\times (1.232444718144824\times {10}^{18}{t}^2\\ +1.298549367102396\times {10}^{21}t-1.748561849689992\times {10}^{22})\\ +5\times {10}^{-19}{e}^{-0.83t}(5.883531005499098\times {10}^{15}{t}^2\\ +3.812288952012609\times {10}^{17}t+5.670981666287240\times {10}^{18})\\ -1.548362123191561\times {10}^{-21}{e}^{-2.92t}(7.675137863758596\times {10}^{16}t\\ -3.260246630067594\times {10}^{17})+5.608298265951419\times {10}^{-26}\\ \times {\mathrm{e}}^{-2.52t}(1.062808766767887\times {10}^{20}{t}^2-1.070860668252633\times {10}^{22}t\\ -3.195993137216618\times {10}^{23})+1.348551711300533\times {10}^{-26}\\ \times {\mathrm{e}}^{-1.26t}(1.277207436992847\times {10}^{20}{t}^3+5.483102259736971\times {10}^{21}{t}^2\\ -3.119452387482340\times {10}^{24}t-1.125583065075483\times {10}^{26})\\ +1.325927750611222\times {10}^{-23}{\mathrm{e}}^{-2.12t}(4.325492835496448\times {10}^{18}{t}^2\\ +1.245072612909103\times {10}^{20}t-5.749959513369393\times {10}^{20}),\end{array}$$

(102)

$$\begin{array}{c}z\left(t\right)={\displaystyle \sum _{i=0}^3}{z}_i\left(t\right)\\ =0.9+2.098254299024467{\mathrm{e}}^{-0.4t}+0.4{\mathrm{e}}^{-0.415t}\sin \mathrm{h}\left(0.015t\right)\\ +0.01953488372093023{\mathrm{e}}^{-0.83t}-0.002666666666666667{\mathrm{e}}^{-0.43t}\\ \times \left(305+6t\right)+0.0009492855151178018{\mathrm{e}}^{-1.23t}\\ -0.001129388038067028{\mathrm{e}}^{-1.66t}-5.897086477312794\times {10}^{-51}{\mathrm{e}}^{-0.83t}\\ \times \left(3.551083966805075\times {10}^{46}t+5.211845325099965\times {10}^{48}\right)\\ +3.944516546299514\times {10}^{-50}{\mathrm{e}}^{-1.26t}(3.552306787033270\times {10}^{45}t\\ -1.231131263500689\times {10}^{46})+1.651476380411790\times {10}^{-48}{\mathrm{e}}^{-0.86t}\\ \times \left(2.162937792347110\times {10}^{46}+6.269476779331703\times {10}^{44}t\right)\\ -1.552529119388551\times {10}^{-61}{\mathrm{e}}^{-0.43t}(4.122305932992456\times {10}^{57}{t}^2\\ +4.191011031876038\times {10}^{59}t+1.422695266527552\times {10}^{61})\end{array}$$

(103)

4.3. Comparison of the Solutions for SIR Epidemic Model

The transmission dynamics of the SIR disease model has been studied numerically. We have coded the New Iterative Method (NIM) and Modified New Iterative Method (MNIM) algorithms for SIR model to test the validity of proposed methods. We have done comparison of absolute errors between ADM [3], NIM and MNIM with respect to RK4 method in

Table 1. Comparisons of absolute errors between ADM, NIM and MNIM solutions with RK4 for SIR Epidemic Model with respect to time $t$ (in years).

$$∆=\left|RK{4}_{h=0.001}-ADM\right|$$				$$∆=\left|RK{4}_{h=0.001}-NIM\right|$$			$$∆=\left|RK{4}_{h=0.001}-MNIM\right|$$
$$\mathit{t}$$	$$∆\mathit{x}\left(\mathit{t}\right)$$	$$∆\mathit{y}\left(\mathit{t}\right)$$	$$∆\mathit{z}\left(\mathit{t}\right)$$	$$∆\mathit{x}\left(\mathit{t}\right)$$	$$∆\mathit{y}\left(\mathit{t}\right)$$	$$∆\mathit{z}\left(\mathit{t}\right)$$	$$∆\mathit{x}\left(\mathit{t}\right)$$	$$∆\mathit{y}\left(\mathit{t}\right)$$	$$∆\mathit{z}\left(\mathit{t}\right)$$
1	0.0001554	0.0001462	0.00853	0.3168	0.01489	0.302	1.119 × 10−4	1.111 × 10−4	1.286 × 10−5
2	0.007425	0.006855	0.1363	0.1776	0.02462	0.2022	3.705 × 10−4	3.637 × 10−4	7.565 × 10−5
3	0.06594	0.05978	0.6892	0.09531	0.03912	0.1348	4.592 × 10−4	4.442 × 10−4	1.518 × 10−4
4	0.3071	0.2745	2.175	0.04445	0.0421	0.0907	4.396 × 10−4	4.188 × 10−4	2.064 × 10−4
5	1.031	0.9136	5.305	0.007361	0.05112	0.06706	4.055 × 10−4	3.817 × 10−4	2.345 × 10−4
6	2.843	2.514	11	0.08268	0.06057	0.06678	3.849 × 10−4	3.598 × 10−4	2.427 × 10−4
7	6.865	6.081	20.37	0.2175	0.3697	0.1133	3.76 × 10−4	3.504 × 10−4	2.379 × 10−4
8	15.02	13.36	34.79	363.9	363.9	0.6809	3.704 × 10−4	3.447 × 10−4	2.249 × 10−4
9	30.38	27.19	55.82	4.648 × 104	4.648 × 104	6.834	3.615 × 10−4	3.358 × 10−4	2.071 × 10−4
10	57.64	51.93	85.3	2.629 × 106	2.629 × 106	54.19	3.459 × 10−4	3.206 × 10−4	1.864 × 10−4

. The MNIM solutions match with RK4 at least 4 decimal places, while both the ADM and NIM solutions are diverging in nature as the value of time $t$ increases. In Figure 2, it can be clear that the 5-iterations of NIM solutions are getting inaccurate as time $t$ gets longer, in contrast the 3-iterations of modified NIM provide excellent accuracy compared to RK4. We have also displayed comparison between RK4 and MNIM on the time interval $\left[0,10\right]$ in Figure 3. The comparisons clearly emphasized the accuracy and validity of the MNIM and RK4, as both outcomes clearly overlap each other over a ten-year time period.

Figure 2a shows that the proportion of susceptible individuals decreases exponentially in the total population. In addition, Figure 2b indicates that the infected population increases suddenly when the epidemic is at its pick level, and after a certain period of time it slow down. Consequently, Figure 2c illustrates that the number of recovered individuals increases rapidly within five years of the epidemic. The impact of high vaccination coverage on the initial population groups with low levels of infection rates, as depicted in Figure 3. The populations of the susceptible and infectious groups decline over time, meanwhile the population of the recovered group grows due to the inclusion of vaccinated and recovered people with permanent immunity, and the disease outbreak stops.

In addition, traditional numerical approaches, such as the Runge–Kutta method, can produce approximate-numeric results. However, analytical series solutions to differential equations can be required. Numerical methods, by nature, cannot yield explicit series solutions. Moreover, the Runge–Kutta method cannot provide the discretized numeric solution in a particular point as it entirely depends on a fixed step size. In contrast the semi-analytical method such as NIM and MNIM can come up with explicit series solutions as well as numerical solutions.

5. Conclusions

This study presents a primer for analyzing and simulating a mathematical model for understanding the dynamics of SIR epidemic model via iterative procedures. Although New Iterative Method (NIM) is an effective method for solving differential equations, this method sometimes fails to handle nonlinear terms from the differential equations. These difficulties may be overcome by coupling the Laplace transform with that of NIM. This combined method offers a rapidly convergent series of solutions. The main characteristics features of Modified New Iterative Method (MNIM) is that it is very simple procedure that does not require any special polynomials or multipliers to be calculated. In particular, it is highly remarkable that the improvement of NIM scheme, i.e., modified NIM provides an excellent accuracy in a large time domain, and it also reduces the computation steps. The effectiveness of MNIM in this research motivates us to extend its applicability to more complex infectious disease models. Further work should reconsider applying these techniques in developing a dynamic model for other conditions such as heart disease [35,36,37], suicide prevention [38], and combine the suggested model with machine learning techniques in developing optimal solutions for infectious diseases such as COVID-19, Pneumonia, and so on [39,40,41].