**1. Introduction**

18 Will-be-set-by-IN-TECH

[4] Esquivel-Flores, O., Benítez-Pérez, H., Méndez Monroy, E. & Menéndez, A. [2010]. Efficient overloading techniques for primary-backup scheduling in real-time systems,

[5] Gupta, R. & Chow, M.-Y. [2010]. Networked control system: Overview and research

[6] Henriksson, D., Cervin, A., Andersson, M. & Årzén, K.-E. [2006]. Truetime: Simulation of networked computer control systems, *Proceedings of the 2nd IFAC Conference on Analysis*

[7] Henriksson, D., Redell, O., El-Khoury, J., Cervin, A., Törngren, M. & Årzén, K.-E. [2006]. Tools for real-time control systems co-design, *in* H. Hansson (ed.), *ARTES – A network for Real-Time research and graduate Education in Sweden 1997–2006*, Department of Information

[8] Lian, F.-L., Moyne, J. & Tilbury, D. [2001]. Time delay modeling and sample time selection for networked control systems. *International Mechanical Engineering Congress and*

[9] Lian, F.-L., Moyne, J. & Tilbury, D. [2002]. Network design consideration for distributed control systems, *Control Systems Technology, IEEE Transactions on* 10(2): 297 –307. [10] Lian, F.-L., Yook, J., Otanez, P., Tilbury, D. & Moyne, J. [2003]. Design of sampling and transmission rates for achieving control and communication performance in networked agent systems, *American Control Conference 2003, Proceedings of* , Vol. 4, pp. 3329 – 3334. [11] Liu, C. L. & Layland, J. W. [1973]. Scheduling algorithms for multiprogramming in a

[12] Ohlin, M., Henriksson, D. & Cervin, A. [2007]. *TRUETIME 1.5 Reference Manual*.

[13] Quanser [2006]. *Quanser 2 DOF Helicopter. User and Control Manual*, Quanser.

[14] Tipsuwan, Y. & Chow, M.-Y. [2003]. Control methodologies in networked control

[15] Xia, F. & Sun, Y. [2008]. *Control and Scheduling Codesign. Flexible Resource Management in*

trends, *Industrial Electronics, IEEE Transactions on* 57(7): 2527 –2535.

*ICI Express Letters Part B: Applications* 1(1): 93–98.

*and Design of Hybrid Systems*, Alghero, Italy.

Technology, Uppsala University, Sweden.

*Exposition, Proceedings of ASME-DSC*, Vol. XX.

hard-real-time environment, *J. ACM* 20(1): 46–61.

systems, *Control Engineering Practice* 11(10): 1099 – 1111.

Innovate-Educate.

*Real-Time Control Systems*, Springer.

Most problems arising from mathematical epidemiology are often described in terms of differential equations. However, it is often very difficult to obtain closed form solutions of such equations, especially those that are nonlinear. In most cases, attempts are made to obtain only approximate or numerical solutions. In this work, we revisit the SIR epidemic model with constant vaccination strategy that was considered in [11], where the Adomian decomposition method was used to solve the governing system of nonlinear initial value differential equations.

In this work we develop new accurate iterative schemes which are based on extending Taylor series based linearization method to obtain accurate and fast converging sequence of hybrid iteration schemes. At first order, the hybrid iteration scheme reduces to quasilinearization method (QLM) which was originally developed in [1]. More recently Mandelzweig and his co-workers [8–10] have extended the application of the QLM to a wide variety of nonlinear BVPs and established that the method converges quadratically. In this work we demonstrate that the proposed hybrid iteration schemes are more accurate and converge faster than the QLM approach.

To implement the method we consider the SIR model that describes the temporal dynamics of a childhood disease in the presence of a preventive vaccine. In SIR models the population is assumed to be divided into the standard three classes namely, the susceptibles (*S*), who can catch the infection but are so far uninfected, the infectives (*I*), those who have the disease and can transmit it to the susceptibles, and the removed (*R*), who have either died or who have recovered and are therefore immune.

The governing equations for the problem are described [11] by

$$\frac{d\mathcal{S}}{dt} = (1 - P)\pi N - \beta \frac{SI}{N} - \mu \mathcal{S}\_{\prime} \tag{1}$$

$$\frac{dI}{dt} = \beta \frac{SI}{N} - (\kappa + \mu)I,\tag{2}$$

©2012 Motsa and Shateyi, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 Will-be-set-by-IN-TECH 68 Numerical Simulation – From Theory to Industry On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology <sup>3</sup>

$$\frac{d\mathcal{R}}{dt} = P\pi N + \kappa I - \mu \mathcal{R}\_\prime \tag{3}$$

*gj*(*z*1, *z*2) = *fj*(*z*1, *z*2) − *fj*(*z*1,*γ*, *z*2,*γ*) −

<sup>1</sup> + *πz*1, *L*2*z*<sup>2</sup> = *z*�

where

where

*Ljzj*,*i*+<sup>1</sup> +

subject to

*L*1*z*<sup>1</sup> = *z*�

We write equation (9) as

differential equation

2 ∑ *s*=1

which can be written as

*zs*,*i*+<sup>1</sup>

which yields the iteration scheme

*∂ fj ∂zs*

*Ljzj*,*i*+<sup>1</sup> +

2 ∑ *s*=1 *zs*,*<sup>i</sup> ∂gj ∂zs*

*Ljzj*,0 +

*Ljzj*,*r*+<sup>1</sup> +

2 ∑ *s*=1 *zs*,0 *∂ fj ∂zs*

2 ∑ *s*=1

*zs*,*r*+<sup>1</sup>

We note that equation (19) is the standard QLM iteration scheme for solving (4 - 5).

of nonlinear equations of the form *f*(*x*) = 0.

*Ljzj* +

2 ∑ *s*=1 *zs ∂ fj ∂zs*

Φ*j*(*z*1,*γ*, *z*2,*γ*) = Ψ*<sup>j</sup>* +

2 ∑ *s*=1

*f*1(*z*1, *z*2) = *βz*1*z*2, *f*2(*z*1, *z*2) = −*βz*1*z*2, Ψ<sup>1</sup> = (1 − *P*)*π*, Ψ<sup>2</sup> = 0. (12)

This idea of introducing the coupled equations of the form (9-10) have previously been used in [3] the construction of Newton-like iteration formulae for the computation of the solutions

We use the quasilinearization method (QLM) of Bellman and Kalaba [1] to solve equation (13). The QLM determines the (*i* + 1)th iterative approximation *zj*,*i*+<sup>1</sup> as the solution of the

(*z*1,*γ*, *z*2,*γ*) +

We assume that *zj*,0 is obtained as a solution of the linear part of equation (13) given by

*∂ fj ∂zs*

2 ∑ *s*=1

> *∂gj ∂zs*

(*z*1,*i*, *z*2,*i*) − *gj*(*z*1,*i*, *z*2,*i*) + Φ*j*(*z*1,*γ*, *z*2,*γ*),

(*zs*,*i*+<sup>1</sup> − *zs*,*i*)

(*z*1,*i*, *z*2,*i*)

2 ∑ *s*=1 *zs*,*<sup>γ</sup> ∂ fj ∂zs*

(*z*1,*γ*, *z*2,*γ*) + *gj*(*z*1,*i*, *z*2,*i*) +

 *∂ fj ∂zs*

2 ∑ *s*=1 (*zs* − *zs*,*γ*)

On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology 69

*∂ fj ∂zs*

<sup>2</sup> + (*π* + *κ*)*z*2, (11)

(*z*1,*γ*, *z*2,*γ*) + *gj*(*z*1, *z*2) = Φ*j*(*z*1,*γ*, *z*2,*γ*), (13)

*∂gj ∂zs*

*z*1,*i*+<sup>1</sup> = *s*0, *z*2,*i*+<sup>1</sup> = *i*0. (17)

(*z*1,*γ*, *z*2,*γ*) = Φ*j*(*z*1,*γ*, *z*2,*γ*), (18)

(*z*1,*r*, *z*2,*r*) = Φ*j*(*z*1,*r*, *z*2,*r*). (19)

(*z*1,*γ*, *z*2,*γ*) − *fj*(*z*1,*γ*, *z*2,*γ*). (14)

(*z*1,*i*, *z*2,*i*) = Φ*j*(*z*1,*γ*, *z*2,*γ*),

*zs*,*i*+<sup>1</sup> = (16)

(15)

(*z*1,*γ*, *z*2,*γ*), (10)

where *S*(*t*), *I*(*t*) and *R*(*t*) denote the susceptibles, infectives and the removed classes respectively.

The total population is denoted by *N* = *S* + *I* + *R*, *μ* is the death rate, *P* is the fraction of citizens vaccinated at birth each year, *β* is the average contact rate, *π* is the constant birth rate, and *κ* is the rate at which an individual recovers from the disease and enters the removed group which also contains vaccinated individuals. Equations (1 - 3) are solved using the new hybrid iteration schemes and the results are compared with results from the Runge-Kutta MATLAB in-built solver ode45.

#### **2. Numerical solution**

To simplify the formulation of the solution, equations (1) - (3) are scaled by dividing by *N*. We define new variables *z*<sup>1</sup> = *S*/*N*, *z*<sup>2</sup> = *I*/*N* and *z*<sup>3</sup> = *R*/*N*. This leads to *z*<sup>1</sup> + *z*<sup>2</sup> + *z*<sup>3</sup> = 1 and if we assume that *π* = *μ*, the scaled new system becomes

$$z\_1'(t) = (1 - P)\pi - \beta z\_1(t)z\_2(t) - \pi z\_1(t), \quad z\_1(0) = s\_{0\prime} \tag{4}$$

$$z\_2'(t) = \beta z\_1(t)z\_2(t) - (\pi + \kappa)z\_2(t), \qquad \qquad z\_2(0) = i\_{0\prime} \tag{5}$$

where *s*<sup>0</sup> and *i*<sup>0</sup> are given constants. The solution for *z*3(*t*) can be obtained from *z*<sup>3</sup> = 1 − *z*<sup>1</sup> − *z*2.

Previous studies [4–7, 12] have shown that the long term behaviour of systems like (4) - (5) can be classified into two categories namely, endemic or eradication. From the long term behaviour of *z*1(*t*) and *z*2(*t*) it holds that the solution asymptotically approaches a disease free equilibrium (DFE) or the endemic equilibrium (EE) where

$$\lim\_{t \to +\infty} (z\_1(t), z\_2(t)) = \text{DFE} = (1 - P\_\prime 0), \tag{6}$$

$$\lim\_{t \to +\infty} (z\_1(t), z\_2(t)) = \text{EE} = \left(\frac{1 - P}{R\_v}, \frac{\pi}{\beta}(R\_v - 1)\right). \tag{7}$$

Here *Rv*, the vaccination reproduction number, is the threshold that determines the stability of the equilibria and is defined by

$$R\_{\mathcal{V}} = \frac{\beta(1 - P)}{\gamma + \pi}. \tag{8}$$

It was shown in [11] that the DFE is locally stable if *Rv <* 1 and the EE is locally stable provided 1 *< Rv* ≤ 4(*κ* + *π*)/*π*. In this work, we use develop new iteration schemes to solve the system (4) - (5) using parameters that yield both the DFE and EE.

#### **3. Method of solution**

To develop the method of solution, we assume that the true solution of (4 - 5) is *zs*,*<sup>α</sup>* (*s* = 1, 2) and *zs*,*<sup>γ</sup>* are the initial approximations. We introduce the following coupled system,

$$(L\_{\dot{\jmath}}z\_{\dot{\jmath}} + f\_{\dot{\jmath}}(z\_{1,\gamma\_{\prime}}z\_{2,\gamma\_{\prime}}) + \sum\_{s=1}^{2} (z\_s - z\_{s,\gamma\_{\prime}}) \frac{\partial f\_{\dot{\jmath}}}{\partial z\_s}(z\_{1,\gamma\_{\prime}}z\_{2,\gamma\_{\prime}}) + g\_{\dot{\jmath}}(z\_1, z\_2) = \Psi\_{\dot{\jmath}} \tag{9}$$

68 Numerical Simulation – From Theory to Industry On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology <sup>3</sup> On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology 69

$$\log\_{\hat{f}}(z\_1, z\_2) = f\_{\hat{f}}(z\_1, z\_2) - f\_{\hat{f}}(z\_{1, \gamma}, z\_{2, \gamma}) - \sum\_{s=1}^{2} (z\_s - z\_{s, \gamma}) \frac{\partial f\_{\hat{f}}}{\partial z\_s}(z\_{1, \gamma\_{\Gamma'}} z\_{2, \gamma}),\tag{10}$$

where

2 Will-be-set-by-IN-TECH

where *S*(*t*), *I*(*t*) and *R*(*t*) denote the susceptibles, infectives and the removed classes

The total population is denoted by *N* = *S* + *I* + *R*, *μ* is the death rate, *P* is the fraction of citizens vaccinated at birth each year, *β* is the average contact rate, *π* is the constant birth rate, and *κ* is the rate at which an individual recovers from the disease and enters the removed group which also contains vaccinated individuals. Equations (1 - 3) are solved using the new hybrid iteration schemes and the results are compared with results from the Runge-Kutta

To simplify the formulation of the solution, equations (1) - (3) are scaled by dividing by *N*. We define new variables *z*<sup>1</sup> = *S*/*N*, *z*<sup>2</sup> = *I*/*N* and *z*<sup>3</sup> = *R*/*N*. This leads to *z*<sup>1</sup> + *z*<sup>2</sup> + *z*<sup>3</sup> = 1 and

where *s*<sup>0</sup> and *i*<sup>0</sup> are given constants. The solution for *z*3(*t*) can be obtained from *z*<sup>3</sup> = 1 − *z*<sup>1</sup> −

Previous studies [4–7, 12] have shown that the long term behaviour of systems like (4) - (5) can be classified into two categories namely, endemic or eradication. From the long term behaviour of *z*1(*t*) and *z*2(*t*) it holds that the solution asymptotically approaches a disease

Here *Rv*, the vaccination reproduction number, is the threshold that determines the stability

*Rv* <sup>=</sup> *<sup>β</sup>*(<sup>1</sup> <sup>−</sup> *<sup>P</sup>*)

It was shown in [11] that the DFE is locally stable if *Rv <* 1 and the EE is locally stable provided 1 *< Rv* ≤ 4(*κ* + *π*)/*π*. In this work, we use develop new iteration schemes to solve

To develop the method of solution, we assume that the true solution of (4 - 5) is *zs*,*<sup>α</sup>* (*s* = 1, 2)

*∂ fj ∂zs*

and *zs*,*<sup>γ</sup>* are the initial approximations. We introduce the following coupled system,

(*zs* − *zs*,*γ*)

<sup>1</sup>(*t*)=(1 − *P*)*π* − *βz*1(*t*)*z*2(*t*) − *πz*1(*t*), *z*1(0) = *s*0, (4)

<sup>2</sup>(*t*) = *βz*1(*t*)*z*2(*t*) − (*π* + *κ*)*z*2(*t*), *z*2(0) = *i*0, (5)

(*z*1(*t*), *z*2(*t*)) = DFE = (1 − *P*, 0), (6)

*<sup>β</sup>* (*Rv* <sup>−</sup> <sup>1</sup>)

*<sup>γ</sup>* <sup>+</sup> *<sup>π</sup>* . (8)

(*z*1,*γ*, *z*2,*γ*) + *gj*(*z*1, *z*2) = Ψ*j*, (9)

. (7)

<sup>1</sup> <sup>−</sup> *<sup>P</sup> Rv* , *π*

*dt* <sup>=</sup> *<sup>P</sup>π<sup>N</sup>* <sup>+</sup> *<sup>κ</sup><sup>I</sup>* <sup>−</sup> *<sup>μ</sup>R*, (3)

*dR*

respectively.

*z*2.

MATLAB in-built solver ode45.

*z*�

*z*�

if we assume that *π* = *μ*, the scaled new system becomes

free equilibrium (DFE) or the endemic equilibrium (EE) where

(*z*1(*t*), *z*2(*t*)) = EE =

the system (4) - (5) using parameters that yield both the DFE and EE.

2 ∑ *s*=1

lim *t*→+∞

lim *t*→+∞

*Ljzj* + *fj*(*z*1,*γ*, *z*2,*γ*) +

of the equilibria and is defined by

**3. Method of solution**

**2. Numerical solution**

$$L\_1 z\_1 = z\_1' + \pi z\_1, \quad L\_2 z\_2 = z\_2' + (\pi + \kappa) z\_2. \tag{11}$$

$$f\_1(z\_1, z\_2) = \beta z\_1 z\_2 \quad f\_2(z\_1, z\_2) = -\beta z\_1 z\_2 \quad \Psi\_1 = (1 - P)\pi \llcorner \Psi\_2 = 0. \tag{12}$$

This idea of introducing the coupled equations of the form (9-10) have previously been used in [3] the construction of Newton-like iteration formulae for the computation of the solutions of nonlinear equations of the form *f*(*x*) = 0.

We write equation (9) as

$$L\_{\hat{f}}z\_{\hat{f}} + \sum\_{s=1}^{2} z\_{s} \frac{\partial f\_{\hat{j}}}{\partial z\_{s}} (z\_{1,\gamma\_{1}}, z\_{2,\gamma\_{1}}) + g\_{\hat{j}}(z\_{1}, z\_{2}) = \Phi\_{\hat{j}}(z\_{1,\gamma\_{1}}, z\_{2,\gamma\_{1}}),\tag{13}$$

where

$$\Phi\_{\dot{\jmath}}(z\_{1,\gamma\prime}z\_{2,\gamma}) = \Psi\_{\dot{\jmath}} + \sum\_{s=1}^{2} z\_{s,\gamma} \frac{\partial f\_{\dot{\jmath}}}{\partial z\_s} (z\_{1,\gamma\prime}z\_{2,\gamma}) - f\_{\dot{\jmath}}(z\_{1,\gamma\prime}z\_{2,\gamma}).\tag{14}$$

We use the quasilinearization method (QLM) of Bellman and Kalaba [1] to solve equation (13). The QLM determines the (*i* + 1)th iterative approximation *zj*,*i*+<sup>1</sup> as the solution of the differential equation

$$L\_j z\_{j, i+1} + \sum\_{s=1}^2 z\_{s, i+1} \frac{\partial f\_j}{\partial z\_s}(z\_{1, \gamma}, z\_{2, \gamma}) + g\_j(z\_{1, i}, z\_{2, i}) + \sum\_{s=1}^2 (z\_{s, i+1} - z\_{s, i}) \frac{\partial g\_j}{\partial z\_s}(z\_{1, \nu}, z\_{2, i}) = \Phi\_j(z\_{1, \gamma}, z\_{2, \gamma}), \tag{15}$$

which can be written as

$$\begin{split} L\_{j}z\_{j,i+1} + \sum\_{s=1}^{2} \left[ \frac{\partial f\_{j}}{\partial z\_{s}}(z\_{1,\gamma\_{i}}z\_{2,\gamma\_{i}}) + \frac{\partial g\_{j}}{\partial z\_{s}}(z\_{1,i},z\_{2,i}) \right] z\_{s,i+1} &= \\ \sum\_{s=1}^{2} z\_{s,i} \frac{\partial g\_{j}}{\partial z\_{s}}(z\_{1,i},z\_{2,i}) - g\_{j}(z\_{1,i},z\_{2,i}) + \Phi\_{j}(z\_{1,\gamma\_{i}},z\_{2,\gamma\_{i}}), \end{split} \tag{16}$$

subject to

$$z\_{1,i+1} = s\_{0\prime} \quad \quad z\_{2,i+1} = i\_0. \tag{17}$$

We assume that *zj*,0 is obtained as a solution of the linear part of equation (13) given by

$$L\_{\vec{\gamma}} z\_{\vec{\gamma},0} + \sum\_{s=1}^{2} z\_{s,0} \frac{\partial f\_{\vec{\gamma}}}{\partial z\_s} (z\_{1,\gamma\_{\prime}}, z\_{2,\gamma\_{\prime}}) = \Phi\_{\vec{\gamma}} (z\_{1,\gamma\_{\prime}}, z\_{2,\gamma\_{\prime}})\_{\prime} \tag{18}$$

which yields the iteration scheme

$$L\_j z\_{j,r+1} + \sum\_{s=1}^2 z\_{s,r+1} \frac{\partial f\_j}{\partial z\_s}(z\_{1,r}, z\_{2,r}) = \Phi\_j(z\_{1,r}, z\_{2,r}).\tag{19}$$

We note that equation (19) is the standard QLM iteration scheme for solving (4 - 5).

#### 4 Will-be-set-by-IN-TECH 70 Numerical Simulation – From Theory to Industry On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology <sup>5</sup>

When *i* = 0 in (16) we can approximate *zj* as

$$z\_{\mathbf{j}} \approx z\_{\mathbf{j},\mathbf{1}}.\tag{20}$$

defined to the interval [-1,1] on which the spectral method can be implemented. We use the transformation *t* = *tF*(*τ* + 1)/2 to map the interval [0, *tF*] to [-1,1], where *tF* is a finite time. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix *D* which is used to approximate the derivatives of the unknown variables *z* at the

where *N* + 1 is the number of collocation points (grid points), **D** = 2*D*/*tF*, and **Z** =

**<sup>a</sup>***ji* <sup>=</sup> *<sup>∂</sup> fi ∂zi*

In this section we present the results of solving the governing equations (4-5) using the iteration scheme-m. For illustration purposes we present the results for *m* = 0, 1, 2 to illustrate the effect of increasing *m* in the accuracy and convergence of the iteration schemes. The number of collocations points in all the results presented here is *N* = 40. In order to assess the accuracy of the proposed method, the present numerical results were compared against results generated using the MATLAB initial value solver ode45. In the numerical simulations presented here, following [11], the governing parameters were carefully selected in order to represent the cases which give rise to both the disease free equilibrium (DFE) and endemic

In this case we observe that *Rv* = 0.186 *<* 1, hence we expect the disease to be eradicated

In this case we observe that *Rv* = 0.186 *<* 1 and as expected, using these parameters, the

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

*z*1,*r*+1(*τ*0) *z*1,*r*+1(*τ*1) . . . *z*1,*r*+1(*τN*) *z*2,*r*+1(*τ*0) *z*2,*r*+1(*τ*1) . . . *z*2,*r*+1(*τN*)

On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology 71

**D***jkz*(*τk*) = **DZ**, *j* = 0, 1, . . . , *N*, (27)

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

Φ1,*r*+1(*τ*0) Φ1,*r*+1(*τ*1) . . . Φ1,*r*+1(*τN*) Φ2,*r*+1(*τ*0) Φ2,*r*+1(*τ*1) . . . Φ2,*r*+1(*τN*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

, (28)

(29)

=

collocation points as the matrix vector product

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

where

*dz dt* <sup>=</sup>

**D** + *πI* + **a**<sup>11</sup> **a**<sup>12</sup>

and *I* is an (*N* + 1) × (*N* + 1) identity matrix.

equilibrium (EE). We consider the following cases

from the population after some time.

disease should be eradicated as *t* → ∞.

1. Case 1: *s*<sup>0</sup> = 1, *i*<sup>0</sup> = 0, *β* = 0.8, *κ* = 0.03, *π* = 0.4, *P* = 0.9.

2. Case 2: *s*<sup>0</sup> = 0.8, *i*<sup>0</sup> = 0.2 *β* = 0.8, *κ* = 0.03, *π* = 0.4, *P* = 0.9.

**4. Results and discussion**

*N* ∑ *k*=0

[*z*(*τ*0), *z*(*τ*1),..., *z*(*τN*)]*<sup>T</sup>* is the vector function at the collocation points *τj*.

Applying the Chebyshev spectral method to (19), for instance, gives

**a**<sup>21</sup> **D** + (*π* + *κ*)*I* + **a**<sup>22</sup>

Thus, setting *i* = 0 in (16) we obtain

$$\begin{split} L\_{j}z\_{j,1} + \sum\_{s=1}^{2} \left[ \frac{\partial f\_{j}}{\partial z\_{s}} (z\_{1,\gamma\_{\prime}}, z\_{2,\gamma}) + \frac{\partial g\_{j}}{\partial z\_{s}} (z\_{1,0}, z\_{2,0}) \right] z\_{s,1} &= \sum\_{s=1}^{2} z\_{s,0} \frac{\partial g\_{j}}{\partial z\_{s}} (z\_{1,0}, z\_{2,0}) \\ & - g\_{j} (z\_{1,0}, z\_{2,0}) + \Phi\_{j} (z\_{1,\gamma\_{\prime}} z\_{2,\gamma}), \end{split} \tag{21}$$

which yields the iteration scheme

$$\begin{split} \mathcal{L}\_{j}\boldsymbol{z}\_{j,r+1} + \sum\_{s=1}^{2} \left[ \frac{\partial f\_{j}}{\partial \boldsymbol{z}\_{s}} (\boldsymbol{z}\_{1,r}, \boldsymbol{z}\_{2,r}) + \frac{\partial g\_{j}}{\partial \boldsymbol{z}\_{s}} (\boldsymbol{z}\_{1,r+1}^{(0)}, \boldsymbol{z}\_{2,r+1}^{(0)}) \right] \boldsymbol{z}\_{s,r+1} &= \sum\_{s=1}^{2} \boldsymbol{z}\_{s,r+1}^{(0)} \frac{\partial g\_{j}}{\partial \boldsymbol{z}\_{s}} (\boldsymbol{z}\_{1,r+1}^{(0)}, \boldsymbol{z}\_{2,r+1}^{(0)}) \\ & \quad - \operatorname{g}\_{j} (\boldsymbol{z}\_{1,r+1}^{(0)}, \boldsymbol{z}\_{2,r+1}^{(0)}) + \boldsymbol{\Phi}\_{j} (\boldsymbol{z}\_{1,r}, \boldsymbol{z}\_{2,r}), \text{ (22)} \end{split} \tag{22}$$

where *z* (0) *<sup>j</sup>*,*r*+<sup>1</sup> is the solution of

$$L\_j z\_{j,r+1}^{(0)} + \sum\_{s=1}^2 z\_{s,r+1}^{(0)} \frac{\partial f\_j}{\partial z\_s}(z\_{1,r}, z\_{2,r}) = \Phi\_j(z\_{1,r}, z\_{2,r}).\tag{23}$$

The general iteration scheme obtained by setting *i* = *m* (*m* ≥ 2) in equation (16), hereinafter referred to as scheme-*m* is

$$\begin{split} \mathcal{L}\_{j}z\_{j,r+1} + \sum\_{s=1}^{2} \left[ \frac{\partial f\_{j}}{\partial z\_{s}}(z\_{1,r}, z\_{2,r}) + \frac{\partial g\_{j}}{\partial z\_{s}}(z\_{1,r+1}^{(m-1)}, z\_{2,r+1}^{(m-1)}) \right] z\_{s,r+1} &= \sum\_{s=1}^{2} z\_{s,r+1}^{(m-1)} \frac{\partial g\_{j}}{\partial z\_{s}}(z\_{1,r+1}^{(m-1)}, z\_{2,r+1}^{(m-1)}) \\ & - g\_{j}(z\_{1,r+1}^{(m-1)}, z\_{2,r+1}^{(m-1)}) + \Phi\_{j}(z\_{1,r}, z\_{2,r}) \end{split} \tag{24}$$

where *z* (*m*−1) *<sup>j</sup>*,*r*+<sup>1</sup> is obtained as the solution of

$$\begin{split} L\_j z\_{j,r+1}^{(m-1)} + \sum\_{s=1}^2 \left[ \frac{\partial f\_j}{\partial z\_s} (z\_{1,r}, z\_{2,r}) + \frac{\partial g\_j}{\partial z\_s} (z\_{1,r+1}^{(m-2)}, z\_{2,r+1}^{(m-2)}) \right] z\_{s,r+1}^{(m-1)} &= \sum\_{s=1}^2 z\_{s,r+1}^{(m-2)} \frac{\partial g\_j}{\partial z\_s} (z\_{1,r+1}^{(m-2)}, z\_{2,r+1}^{(m-2)}) \\ &- g\_j (z\_{1,r+1}^{(m-2)}, z\_{2,r+1}^{(m-2)}) + \Phi\_j (z\_{1,r}, z\_{2,r}). \end{split} \tag{25}$$

The initial approximation for solving the iteration algorithms, scheme-*m* is obtained by solving the linear part of the governing equations (4 - 5). This gives

$$z\_{1,0} = (1 - P)(1 - e^{-\pi t}) + s\_0 e^{-\pi t}, \qquad z\_{2,0} = i\_0 e^{-(\pi + \kappa)t}.\tag{26}$$

The iteration schemes (19),(24 - 25) can be solved numerically using standard methods such as finite difference, finite elements, spline collocation methods,etc. In this study we use the Chebyshev spectral collocation method to solve the iteration schemes. For brevity, we omit the details of the spectral methods, and refer interested readers to ([2, 13]). Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval [-1,1] on which the spectral method can be implemented. We use the transformation *t* = *tF*(*τ* + 1)/2 to map the interval [0, *tF*] to [-1,1], where *tF* is a finite time. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix *D* which is used to approximate the derivatives of the unknown variables *z* at the collocation points as the matrix vector product

$$\frac{dz}{dt} = \sum\_{k=0}^{N} \mathbf{D}\_{jk} z(\tau\_k) = \mathbf{D} \mathbf{Z}, \quad j = 0, 1, \dots, N,\tag{27}$$

where *N* + 1 is the number of collocation points (grid points), **D** = 2*D*/*tF*, and **Z** = [*z*(*τ*0), *z*(*τ*1),..., *z*(*τN*)]*<sup>T</sup>* is the vector function at the collocation points *τj*.

Applying the Chebyshev spectral method to (19), for instance, gives

$$\begin{bmatrix} \mathbf{D} + \boldsymbol{\pi}I + \mathbf{a}\_{11} \\\\ \hline \\ & \\ \hline \\ & \mathbf{a}\_{21} \end{bmatrix} \quad \mathbf{D} + (\boldsymbol{\pi} + \mathbf{\kappa})I + \mathbf{a}\_{22} \\\\ \begin{bmatrix} \mathbf{D} + (\boldsymbol{\pi} + \mathbf{\kappa})I + \mathbf{a}\_{21} \\\\ \hline \\ \mathbf{D} + (\boldsymbol{\pi} + \mathbf{\kappa})I + \mathbf{a}\_{22} \\\\ \hline \\ \mathbf{z}\_{2,r+1}(\tau\_{0}) \\\\ \hline \\ \mathbf{z}\_{2,r+1}(\tau\_{N}) \end{bmatrix} \begin{bmatrix} z\_{1,r+1}(\tau\_{0}) \\\\ \hline \\ z\_{1,r+1}(\tau\_{1}) \\\\ \hline \\ \mathbf{0}\_{2,r+1}(\tau\_{0}) \\\\ \hline \\ \mathbf{0}\_{2,r+1}(\tau\_{1}) \\\\ \hline \\ \mathbf{0}\_{2,r+1}(\tau\_{N}) \end{bmatrix},\tag{28}$$

where

4 Will-be-set-by-IN-TECH

(*z*1,0, *z*2,0)

 *zs*,1 =

*zs*,*r*+<sup>1</sup> =

−*gj*(*z* (0) 1,*r*+1, *z* (0)

*zs*,*r*+<sup>1</sup> =

−*gj*(*z*

2 ∑ *s*=1 *z* (*m*−1) *s*,*r*+1

(*m*−1) 1,*r*+<sup>1</sup> , *z*

2 ∑ *s*=1 *z* (*m*−2) *s*,*r*+1

(*m*−2) 1,*r*+<sup>1</sup> , *z*

−*gj*(*z*

, *z*2,0 = *i*0*e*−(*π*+*κ*)*<sup>t</sup>*

*zj* ≈ *zj*,1. (20)

(*z*1,0, *z*2,0)

−*gj*(*z*1,0, *z*2,0) + Φ*j*(*z*1,*γ*, *z*2,*γ*), (21)

2,*r*+1) + Φ*j*(*z*1,*r*, *z*2,*r*), (22)

2,*r*+<sup>1</sup> ) + Φ*j*(*z*1,*r*, *z*2,*r*),

2,*r*+<sup>1</sup> ) + Φ*j*(*z*1,*r*, *z*2,*r*).

. (26)

(*m*−1) 2,*r*+<sup>1</sup> )

(24)

(25)

(*m*−2) 2,*r*+<sup>1</sup> )

*∂gj ∂zs* (*z* (0) 1,*r*+1, *z* (0) 2,*r*+1)

(*z*1,*r*, *z*2,*r*) = Φ*j*(*z*1,*r*, *z*2,*r*). (23)

*∂gj ∂zs* (*z* (*m*−1) 1,*r*+<sup>1</sup> , *z*

(*m*−1)

*∂gj ∂zs* (*z* (*m*−2) 1,*r*+<sup>1</sup> , *z*

(*m*−2)

2 ∑ *s*=1 *zs*,0 *∂gj ∂zs*

2 ∑ *s*=1 *z* (0) *s*,*r*+1

When *i* = 0 in (16) we can approximate *zj* as

 *∂ fj ∂zs*

(*z*1,*r*, *z*2,*r*) +

*Ljz* (0) *<sup>j</sup>*,*r*+<sup>1</sup> +

(*z*1,*r*, *z*2,*r*) +

*<sup>j</sup>*,*r*+<sup>1</sup> is obtained as the solution of

(*z*1,*r*, *z*2,*r*) +

(*z*1,*γ*, *z*2,*γ*) +

*∂gj ∂zs* (*z* (0) 1,*r*+1, *z* (0) 2,*r*+1) 

2 ∑ *s*=1 *z* (0) *s*,*r*+1

*∂gj ∂zs* (*z* (*m*−1) 1,*r*+<sup>1</sup> , *z*

*∂gj ∂zs* (*z* (*m*−2) 1,*r*+<sup>1</sup> , *z*

solving the linear part of the governing equations (4 - 5). This gives

*<sup>z</sup>*1,0 = (<sup>1</sup> <sup>−</sup> *<sup>P</sup>*)(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*−*π<sup>t</sup>*

*∂ fj ∂zs*

The general iteration scheme obtained by setting *i* = *m* (*m* ≥ 2) in equation (16), hereinafter

(*m*−1) 2,*r*+<sup>1</sup> ) 

(*m*−2) 2,*r*+<sup>1</sup> ) *z* (*m*−1) *<sup>s</sup>*,*r*+<sup>1</sup> =

The initial approximation for solving the iteration algorithms, scheme-*m* is obtained by

) + *s*0*e*−*π<sup>t</sup>*

The iteration schemes (19),(24 - 25) can be solved numerically using standard methods such as finite difference, finite elements, spline collocation methods,etc. In this study we use the Chebyshev spectral collocation method to solve the iteration schemes. For brevity, we omit the details of the spectral methods, and refer interested readers to ([2, 13]). Before applying the spectral method, it is convenient to transform the domain on which the governing equation is

*∂gj ∂zs*

Thus, setting *i* = 0 in (16) we obtain

2 ∑ *s*=1

which yields the iteration scheme

 *∂ fj ∂zs*

*<sup>j</sup>*,*r*+<sup>1</sup> is the solution of

*Ljzj*,1 +

2 ∑ *s*=1

referred to as scheme-*m* is

2 ∑ *s*=1

(*m*−1)

2 ∑ *s*=1  *∂ fj ∂zs*

 *∂ fj ∂zs*

*Ljzj*,*r*+<sup>1</sup> +

(0)

where *z*

*Ljzj*,*r*+<sup>1</sup> +

where *z*

*Ljz* (*m*−1) *<sup>j</sup>*,*r*+<sup>1</sup> +

$$\mathbf{a}\_{ji} = \frac{\partial f\_i}{\partial z\_i} \tag{29}$$

and *I* is an (*N* + 1) × (*N* + 1) identity matrix.

#### **4. Results and discussion**

In this section we present the results of solving the governing equations (4-5) using the iteration scheme-m. For illustration purposes we present the results for *m* = 0, 1, 2 to illustrate the effect of increasing *m* in the accuracy and convergence of the iteration schemes. The number of collocations points in all the results presented here is *N* = 40. In order to assess the accuracy of the proposed method, the present numerical results were compared against results generated using the MATLAB initial value solver ode45. In the numerical simulations presented here, following [11], the governing parameters were carefully selected in order to represent the cases which give rise to both the disease free equilibrium (DFE) and endemic equilibrium (EE). We consider the following cases


#### 6 Will-be-set-by-IN-TECH 72 Numerical Simulation – From Theory to Industry On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology <sup>7</sup>


The results for Case 1 are shown on Figs. 1 - 2. In this case, the initial guess and the first few iterations match the numerical solution all the iterative schemes in the plots of *s*(*t*),*r*(*t*). We observe that *s*(*t*) decreases monotonically with time while *r*(*t*) increases with time. The graph of the profile for *i*(*t*) is not shown because *i*(*t*) = 0 in this particular case.

0 10 20 30

0 10 20 30

Case 1: Scheme−1

On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology 73

ode45 r = 1 r = 2 r = 3 r = 4

t

0 0.2 0.4 0.6 0.8 1

> ode45 r = 1 r = 2 r = 3 r = 4

Case 1: Scheme−2

0 10 20 30

t

0 10 20 30

Case 2: Scheme−2

s(t)

ode45 r = 1 r = 2 r = 3 r = 4

0 0.2 0.4 0.6 0.8

> ode45 r = 1 r = 2 r = 3 r = 4

0 10 20 30

Case 2: Scheme−1

ode45 r = 1 r = 2 r = 3 r = 4

t

t

**Figure 3.** Case 2: Comparison of the numerical solution of the population fractions *s*(*t*) against the

0 0.2 0.4 0.6 0.8

s(t)

0 10 20 30

Case 2: Scheme−0

t

**Figure 2.** Case 1: Comparison of the numerical solution of the population fractions *r*(*t*) against the

r(t)

ode45 r = 1 r = 2 r = 3 r = 4

Case 1: Scheme−0

t

r(t)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

results from the iteration schemes-0, 1 and 2

0 0.2 0.4 0.6 0.8

results from the iteration schemes-0, 1 and 2

s(t)

r(t)

**Figure 1.** Case 1: Comparison of the numerical solution of the population fractions *s*(*t*) against the results from the iteration schemes-0, 1 and 2

Figs. 3 - 5 show the numerical approximation of the profiles of the different classes for Case 2. Again, all the iterative schemes rapidly converge to the numerical solution. The population of the susceptibles decreases with time and that of the removed (those recovered with immunity) increases with time. The infected population initially increases and reaches a maximum, then gradually decreases to zero as *t* → ∞.

Figs. 6 - 8 show the numerical approximation of the profiles of the different classes for Case 3. It can be noted from the graphs that the Scheme-2 converges fastest towards the numerical results. Only 10 iterations are required for full convergence in Scheme-2 compared to 14 iterations in Scheme-1 and 28 iterations in Scheme-1.

Figs. 8 - 11 shows the variation all the population groups with time for Case 4. Again, we observe that Scheme-2 converges fastest towards the numerical results. Only 5 iterations are required for full convergence in Scheme-2 compared to 6 iterations in Scheme-1 and 12 iterations in Scheme-1.

6 Will-be-set-by-IN-TECH

In this case*Rv* = 1.86 *>* 1 which leads to the endemic equilibrium (no disease eradication).

In this case *Rv* = 1.32 *>* 1 which leads to the endemic equilibrium (no disease eradication).

0 0.2 0.4 0.6 0.8

> ode45 r = 1 r = 2 r = 3 r = 4

Case 1: Scheme−2

0 10 20 30

t

Figs. 3 - 5 show the numerical approximation of the profiles of the different classes for Case 2. Again, all the iterative schemes rapidly converge to the numerical solution. The population of the susceptibles decreases with time and that of the removed (those recovered with immunity) increases with time. The infected population initially increases and reaches a maximum, then

Figs. 6 - 8 show the numerical approximation of the profiles of the different classes for Case 3. It can be noted from the graphs that the Scheme-2 converges fastest towards the numerical results. Only 10 iterations are required for full convergence in Scheme-2 compared to 14

Figs. 8 - 11 shows the variation all the population groups with time for Case 4. Again, we observe that Scheme-2 converges fastest towards the numerical results. Only 5 iterations are required for full convergence in Scheme-2 compared to 6 iterations in Scheme-1 and 12

**Figure 1.** Case 1: Comparison of the numerical solution of the population fractions *s*(*t*) against the

0 10 20 30

Case 1: Scheme−1

ode45 r = 1 r = 2 r = 3 r = 4

t

The results for Case 1 are shown on Figs. 1 - 2. In this case, the initial guess and the first few iterations match the numerical solution all the iterative schemes in the plots of *s*(*t*),*r*(*t*). We observe that *s*(*t*) decreases monotonically with time while *r*(*t*) increases with time. The graph

s(t)

ode45 r = 1 r = 2 r = 3 r = 4

3. Case 3: *s*<sup>0</sup> = 0.8, *i*<sup>0</sup> = 0.2 *β* = 0.8, *κ* = 0.03, *π* = 0.4, *P* = 0.

4. Case 4: *s*<sup>0</sup> = 0.8, *i*<sup>0</sup> = 0.2 *β* = 0.8, *κ* = 0.03, *π* = 0.4, *P* = 0.3.

of the profile for *i*(*t*) is not shown because *i*(*t*) = 0 in this particular case.

0 10 20 30

Case 1: Scheme−0

t

s(t)

0 0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8

results from the iteration schemes-0, 1 and 2

gradually decreases to zero as *t* → ∞.

iterations in Scheme-1.

iterations in Scheme-1 and 28 iterations in Scheme-1.

s(t)

**Figure 2.** Case 1: Comparison of the numerical solution of the population fractions *r*(*t*) against the results from the iteration schemes-0, 1 and 2

**Figure 3.** Case 2: Comparison of the numerical solution of the population fractions *s*(*t*) against the results from the iteration schemes-0, 1 and 2

**Figure 4.** Case 2: Comparison of the numerical solution of the population fractions *i*(*t*) against the results from the iteration schemes-0, 1 and 2

0 10 20 30 40 50

Case 3: Scheme−2

0.4

0.6

s(t)

ode45 r = 10 r = 15 r = 20 r = 28

0.8

1

ode45 r = 2 r = 3 r = 4 r = 5

0 10 20 30 40 50

Case 3: Scheme−1

On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology 75

ode45 r = 4 r = 8 r = 12 r = 14

t

t

0 10 20 30 40 50

Case 3: Scheme−2

0

0.2

0.4

i(t)

ode45 r = 2 r = 4 r = 8 r = 10

0 10 20 30 40 50

Case 3: Scheme−1

ode45 r = 4 r =8 r = 12 r = 14

t

t

**Figure 7.** Case 3: Comparison of the numerical solution of the population fractions *i*(*t*) against the

**Figure 6.** Case 3: Comparison of the numerical solution of the population fractions *s*(*t*) against the

ode45 r = 10 r = 15 r = 20 r = 28

0.4

0 10 20 30 40 50

Case 3: Scheme−0

t

i(t)

0

0.2

0.4

0.6

0.6

s(t)

0.8

1

0 10 20 30 40 50

Case 3: Scheme−0

t

0.4

results from the iteration schemes-0, 1 and 2

0

results from the iteration schemes-0, 1 and 2

0.2

0.4

i(t)

0.6

s(t)

0.8

1

**Figure 5.** Case 2: Comparison of the numerical solution of the population fractions *r*(*t*) against the results from the iteration schemes-0, 1 and 2

8 Will-be-set-by-IN-TECH

0 10 20 30

Case 2: Scheme−2

0

0.1

0.2

i(t)

ode45 r = 1 r = 2 r = 3 r = 4

> ode45 r = 1 r = 2 r = 3 r = 4

0 10 20 30

Case 2: Scheme−1

ode45 r = 1 r = 2 r = 3 r = 4

t

t

0 10 20 30

Case 2: Scheme−2

r(t)

0 0.2 0.4 0.6 0.8 1

> ode45 r = 1 r = 2 r = 3 r = 4

0 10 20 30

Case 2: Scheme−1

ode45 r = 1 r = 2 r = 3 r = 4

t

t

**Figure 5.** Case 2: Comparison of the numerical solution of the population fractions *r*(*t*) against the

0 0.2 0.4 0.6 0.8 1

r(t)

0 10 20 30

Case 2: Scheme−0

t

**Figure 4.** Case 2: Comparison of the numerical solution of the population fractions *i*(*t*) against the

ode45 r = 1 r = 2 r = 3 r = 4

0

0.1

0.2

i(t)

0.3

0 10 20 30

Case 2: Scheme−0

t

0

results from the iteration schemes-0, 1 and 2

0 0.2 0.4 0.6 0.8 1

results from the iteration schemes-0, 1 and 2

r(t)

0.1

0.2

i(t)

**Figure 6.** Case 3: Comparison of the numerical solution of the population fractions *s*(*t*) against the results from the iteration schemes-0, 1 and 2

**Figure 7.** Case 3: Comparison of the numerical solution of the population fractions *i*(*t*) against the results from the iteration schemes-0, 1 and 2

**Figure 8.** Case 3: Comparison of the numerical solution of the population fractions *r*(*t*) against the results from the iteration schemes-0, 1 and 2

0 10 20 30 40 50

Case 4: Scheme−2 ode45 r = 2 r = 3 r = 4 r = 5

i(t)

0 0.1 0.2 0.3 0.4 0.5

0 10 20 30 40 50

Case 4: Scheme−1

On New High Order Iterative Schemes for Solving Initial Value Problems in Epidemiology 77

ode45 r = 2 r = 3 r = 4 r = 6

t

t

0 10 20 30 40 50

Case 4: Scheme−2

r(t)

ode45 r = 2 r = 4 r = 6 r = 10

−0.1 0 0.1 0.2 0.3 0.4

> ode45 r = 2 r = 3 r = 4 r = 5

0 10 20 30 40 50

Case 4: Scheme−1

ode45 r = 2 r = 3 r = 4 r = 6

t

t

**Figure 11.** Case 4: Comparison of the numerical solution of the population fractions *r*(*t*) against the

−0.1 0 0.1 0.2 0.3 0.4

r(t)

0 10 20 30 40 50

Case 4: Scheme−0

t

**Figure 10.** Case 4: Comparison of the numerical solution of the population fractions *i*(*t*) against the

0 0.1 0.2 0.3 0.4 0.5

i(t)

0 10 20 30 40 50

Case 4: Scheme−0 ode45 r = 2 r = 6 r = 10 r = 12

t

0

results from the iteration schemes-0, 1 and 2

−0.1 0 0.1 0.2 0.3 0.4

results from the iteration schemes-0, 1 and 2

r(t)

0.2

0.4

i(t)

**Figure 9.** Case 4: Comparison of the numerical solution of the population fractions *s*(*t*) against the results from the iteration schemes-0, 1 and 2

10 Will-be-set-by-IN-TECH

0 10 20 30 40 50

Case 3: Scheme−2

0

0.02

r(t)

0.04

ode45 r = 2 r = 4 r = 8 r = 10

0 10 20 30 40 50

Case 3: Scheme−1

ode45 r = 4 r = 8 r = 12 r = 14

t

0 10 20 30 40 50

Case 4: Scheme−1 ode45 r = 2 r = 3 r = 4 r = 6

t

t

0 10 20 30 40 50

Case 4: Scheme−2 ode45 r = 2 r = 3 r = 4 r = 5

0.4

0.6

s(t)

0.8

1

t

**Figure 9.** Case 4: Comparison of the numerical solution of the population fractions *s*(*t*) against the

**Figure 8.** Case 3: Comparison of the numerical solution of the population fractions *r*(*t*) against the

ode45 r = 10 r = 15 r = 20 r = 28

0 10 20 30 40 50

Case 3: Scheme−0

t

0

0.4

0.6

s(t)

0.8

1

0 10 20 30 40 50

Case 4: Scheme−0 ode45 r = 2 r = 6 r = 10 r = 12

t

0.02

r(t)

0.04

0

results from the iteration schemes-0, 1 and 2

0.4

results from the iteration schemes-0, 1 and 2

0.6

s(t)

0.8

1

0.02

r(t)

0.04

**Figure 10.** Case 4: Comparison of the numerical solution of the population fractions *i*(*t*) against the results from the iteration schemes-0, 1 and 2

**Figure 11.** Case 4: Comparison of the numerical solution of the population fractions *r*(*t*) against the results from the iteration schemes-0, 1 and 2
