EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 4, 2024, 2763-2799 ISSN 1307-5543 – ejpam.com Published by New York Business Global On the Solution of the Fractional-Order Pneumonia Model Using Numerical Computational Methods Ahmad Alalyani Department of Mathematics, Faculty of Science, Al-Baha University, Al-Baha, Saudi Arabia Abstract. This paper deals with the solution and the dynamics of the pneumonia fractional-order mathematical model using numerical computational methods. We study positivity, boundedness, equilibria, local and global stability, and the basic reproductive number R0 of the proposed model, which is the one most significant parameter in epidemiological modeling. It estimates the average number of additional infections induced by a sole infectious individual within a fully susceptible group during the mean period of infection. To verify the theoretical analysis of the proposed model, we use numerical techniques including the Adams-Bashforth-Moulton method, the generalized Euler method, the generalized Runge-Kutta method, and the multistep generalized differential transform method. The numerical results and simulations confirm the convergence between the presented fractional-order model and its integer-order form. The proposed model proves to be a valuable tool for investigating dynamical and numerical analysis for a variety of disease models in epidemiology. 2020 Mathematics Subject Classifications: 26A33, 34A08, 34D20, 65L99, 92-10, 92C42, 92D25, 92D30 Key Words and Phrases: Pneumonia disease, Fractional-order mathematical model, Caputo fractional derivative, Stability analysis, Numerical techniques 1. Introduction According to Johns Hopkins Medicine [34], pneumonia is a disease of one or both of the lungs caused by infection with bacteria, viruses, or fungi. Pneumonia has more than 30 different causes, which are grouped according to the cause. The main species of pneumonia include bacterial pneumonia, viral pneumonia, mycoplasma pneumonia, and other pneumonias. People of different ages can get pneumonia, but the groups at the highest risk are children under the age of 2, adults ages 65 and older, individuals with specific medical conditions, and individuals who smoke. The majority of people with pneumonia have a good reaction to medicine, but pneumonia can be lethal. In 2013, Ong’ala et al. [46] created a mathematical model for pneumonia dynamics with carriers. DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5396 Email address: azaher@bu.edu.sa (A. Alalyani) https://www.ejpam.com 2763 Copyright: © 2024 The Author(s). (CC BY-NC 4.0) A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2764 They use bifurcation and stability of equilibrium points to study the reduction paths or transfer rates between infected individuals and the carriers. In 2014, a simple model of ODE was presented by Mochan et al. [36] to dynamically describe the host’s immune response to pneumonia caused by bacteria infection in murine strains. In 2014, Drusano et al. [21] examined how granulates work to kill bacteria and found no role for antibiotics. Ndelwa et al. [39] created in 2015 a mathematical model to describe the pneumonia transmission process in order to estimate transmission and effects. In 2015, Kosasih et al. [29] used wavelet-based crackle detection to analyze a mathematical model of cough sounds for quick identification of pneumonia caused by bacteria in children. Cesar et al. [16] mathematically evaluated in 2016 a fine particulate matter using a model and assessed medicines for pediatric asthma and pneumonia. Atypical bacterial pathogens were recorded as the primary reasons for lower respiratory illnesses such as CAP, bronchitis, and coughs by Marchello et al. [32] in 2016. Cheng et al. in 2017 estimated an IAV-SP model dynamically and mathematically [17]. To develop the respiratory health of COPD patients, they created a quantitative risk-assessment structure. In 2017, Kosasih and Abeyratne [28] proposed a straightforward mathematical model to illustrate the measurement analysis for the clinical identification of pediatric pneumonia. Tilahun et al. [55] presented in 2017 a non-linear mathematical model along with an analysis of optical control approaches for pneumonia disease. Tilahun et al. [56] consid- ered in 2018 a co-infection model for the diseases of pneumonia and typhoid fever, and their distinctive relation in the case of a treatment and medical plans of action was math- ematically analyzed. Based on mathematical characteristics of cough sounds, Raj et al. in 2018 [49] studied the classification of pneumonia and asthma in poor populations. Kizito and Tumwiine [27] provided a mathematical model in 2018 that illustrates how bacteria control the spread of pneumonia. In addition, the dynamics of vaccine formulation and treatment were discussed. A non-linear mathematical model that describes the modeling of co-infection of the influenza A virus and pneumonia within the host was studied in 2018 by Mbabazi et al. [33]. Tilahun in 2019 [54] used theorems and ordinary differential equations to describe a pneumonia-meningitis co-infection model. Tilahun provided an explanation of different disease clearance methods. In 2019, Diah and Aziz [19] reviewed dynamic mathematical models of pneumonia, followed by earlier research. Tilahun in 2019 [53] investigated a mathematical model of co-dynamics for the diseases meningitis and pneumonia. In 2020, Otoo et al. [47] studied a model that illustrates how bacteria spread pneumonia, and the analysis in this study was helpful in determining the effects of vaccination on disease control. In 2020, the graphical result for the dynamics of a math- ematical model for pneumonia was provided by Zephaniah et al. [60]. Ming et al. [35] discussed in 2020 the growing number of coronavirus pneumonia cases and the spread of this disease in Wuhan, China. In 2020, Jung et al. [26] exhibited the observations through a variety of clinical tests and revealed a novel pathogen as the cause of disease. In 2021, Wafula et al. [59] established a deterministic mathematical model of pneumonia-HIV co- infection with the utilization of anti-pneumonia and ART therapy strategies as control to describe optimal control treatment. Oluwatobi and Erinle-Ibrahim in 2021 [45] inves- tigated the impact of treating pneumonia disease, the existence of efficient reproduction A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2765 numbers, and the stability of equilibrium points. In 2021, a delayed mathematical model- ing technique was implemented by Muhammad Naveed et al. [38] to analyze pneumonia disease. Olumuyiwa James Peter et al. in 2021 [44] introduced a model for pneumococcal pneumonia infection using fractional order derivatives in the sense of the Caputo-Fabrizio operator. In 2022, Sayed Saber et al. [52] studied a fractional mathematical SVCIR model to estimate the pneumonia disease transmission dynamics in sheep and goats in Al-Baha region of Saudi Arabia with cost-effective strategies. The research conducted by Sayed Saber et al. investigated critical parameters that impact the transmission of pneumonia in livestock. They provided evidence that the implementation of cost-effective vaccination and treatment strategies led to a substantial reduction in disease prevalence. Furthermore, the fractional models usually have more accurate predictions compared to the traditional integer models. In 2022, Kamaledin Abodayeh et al. [1] discussed the dynamical and nu- merical analysis of a pneumonia model using appropriate numerical techniques. Fractional epidemic models provide various advantages, such as improved accuracy in representing the memory and hereditary properties of disease transmission. They offer improved ca- pacity to explain the complex dynamics of disease transmission over time, a facet that traditional integer models may neglect. This results in heightened predictive accuracy and more efficacious intervention strategies. There is also substantial literature dealing with the mathematical investigation of a disease model; (see, for example, but not limited to, ([13], [5], [50], [7], [3], [14], [2], [12]). We develop in this research a fractional-order model for the presented model in [1] to describe the dynamics of transmission for pneumonia disease. More precisely, for any arbitrary time t, a fractional-order SCIR model is taken into consideration in the context of the Caputo derivative Dα, 0 < α ≤ 1, and is provided by: DαS(t) = Λ− δ(I(t) + ωC(t)) N S(t)− µS(t) + ηR(t), DαC(t) = δ(I(t) + ωC(t)) N θS(t)− Z1C(t), DαI(t) = δ(I(t) + ωC(t)) N (1− θ)S(t) + πC(t)− Z2I(t), DαR(t) = βC(t) + τI(t)− (µ+ η)R(t), (1) with the initial conditions S(0) = S0, C(0) = C0, I(0) = I0, and R(0) = R0. Here, the infection force is give by δ(I(t)+ωC(t)) N . The fractional-order SCIR model, as expressed in Eq. (1), serves as a mathematical representation tailored to elucidate the dynamics of pneumonia transmission. The model employs the Caputo derivative Dα, where 0 < α ≤ 1, to incorporate memory effects in the system dynamics. This approach facilitates a more precise characterization of biological processes in the real-world, frequently exhibited by fractional dynamics. Fractional mathematical modeling has been shown to be an effective mathematical tool with uses in a variety of fields of study (see, for example, but not limited to, ([4], [8], [10], [11], [48], [9], [51], [57], [6]). In Eq. (1), S(t) represents the susceptible who is susceptible to infection with pneumonia. C(t) stands for the individuals who carry pneumonia bacteria and can transmit the infection. I(t) stands for the infectious A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2766 individuals who have the potential to transfer the infection to the susceptible; and R(t) stands for the individuals who received treatment for pneumonia and made a recovery. Over time, susceptibility to the disease may be regained by these individuals, contingent upon factors such as waning immunity. Λ represents the per capita recruitment rate of new susceptible individuals into the population, which may occur through processes such as births or immigration. δ stands for the transmission rate of pneumonia, denoting the efficiency with which the disease spreads from infectious and carrier individuals to susceptible individuals. ω indicates the rate of treated individuals who had vaccinations. µ stands for the per capita natural death rate of individuals due to non-pneumonia-related causes. η represents the rate of immunity loss in recovered individuals, leading to their return to the susceptible class and accounting for the potential of reinfection. Z1 = µ+β+π is a combined rate representing the sum of (µ), the per capita carriers recovery rate (β), and the rate of developing symptoms by carriers (π). Z2 = τ + µ + σ is a combined rate representing the sum of (µ), the per capita recovery rate of infected individuals by pneumonia (τ), and the disease-induced death rate and birth rate of the human population per capita (σ). θ is the rate of susceptible individuals who, following infection, become carriers instead of manifesting symptoms. The total population size, often assumed to be constant, is denoted by N , where N = S(t) + C(t) + I(t) +R(t). We study the local and global stability of steady states. We derive the fundamental reproduction threshold parameter to show that the disease-free steady state is locally and globally asymptotically stable when R0 < 1. However, a positive (endemic) steady state that is both locally and globally asymptotically stable exists when R0 > 1. To obtain the approximate solution and verify the theoretical analysis of the fractional-order SCIR model, we have used numerical methods including the Adams-Bashforth-Moulton Method (ABMM), the Generalized Euler Method (GEM), the Generalized Runge-Kutta Method (GRKM), and the Multistep Generalized Differential Transform Method (MSGDTM). The paper is structured as follows: Section 2 discusses the properties of the model. Section 3 investigates the stability analysis of the model. Section 4 explains the numerical compu- tational methods of solution. Sections 5, 6, and 7 present the numerical simulations and results, discussion of the results, and conclusion. 2. Properties of the Model Let R+ be the set of all nonnegative real numbers and assume that Ω+ = {(S,C, I,R) ∈ R4 + : S ≥ 0, C ≥ 0, I ≥ 0, R ≥ 0, max (|S|, |C|, |I|, |R|) ≤ N}. Theorem 1 ([13]). If the starting condition Γ0 = (S(0), C(0), I(0), R(0)) ∈ Ω is satisfied, then the solution Γ = (S(t), C(t), I(t), R(t)) ∈ Ω to the model (1) is the only one that can exist for t ≥ 0. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2767 Proof. Take T1(Γ) = Λ− δI(t)S(t) N − δωC(t)S(t) N − µS(t) + δR(t), T2(Γ) = δθI(t)S(t) N + δωθC(t)S(t) N − Z1C(t), T3(Γ) = δ(I(t) + ωC(t)) N (1− θ)S(t) + πC(t)− Z2I(t), T4(Γ) = βC(t) + τI(t)− (µ+ η)R(t). Inputting the values Γ,Γ ∈ Ω, with | δ(I(t)+ωC(t)) N | ≤ K, yields the following: ∥T(Γ)− T(Γ)∥ = |T1(Γ)− T1(Γ)|+ |T2(Γ)− T2(Γ)|+ |T3(Γ)− T3(Γ)|+ |T4(X)− T4(X)| ≤ (3K + µ−Kθ)|S − S|+ (Z1 + β + π)|C − C| + (τ + Z2)|I − I|+ (δ + µ+ η)|R−R| ≤ Γ|X −X|, where Γ = max {(2δ(1 + ω) + β + Z1), (2δω + π + µ), (2δ + τ + Z2), (δ + µ+ η)} . Consequently, the solution to the model (1) exists and is unique iff H(X) satisfies the Lipschitz condition. Lemma 1 ([31], Lemma 3). Let u(t) be a function meeting the conditions that Dαu(t) exists for every t. Therefore, { Dαu(t) ≤ −λu(t) + µ, u(t0) = ut0 , where 0 < α < 1, (λ, µ) ∈ R2, λ ̸= 0, and t0 ≥ 0 is the initial time. Then, u(t) ≤ ( u(t0)− µ λ ) Eα,1 [ − λ(t− t0) α ] + µ λ , where Eα,1(t) = ∞∑ k=0 tk Γ(kq + 1) > 0 is the Mittag Leffler function. Theorem 2 ([31]). The model (1) has non-negative solutions if and only if θ < 1. Proof. We have DαS(t)|S=0 = Λ+ δR(t) > 0, DαC(t)|C=0 = δθI(t)S(t) N + δωθC(t)S(t) N > 0, DαI(t)|I=0 = δΩC(t) N (1− θ)S(t) + πC(t) > 0, DαR(t)|R=0 = βC(t) + τI(t) > 0. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2768 Using Lemmas 5 and 6 from [15], we see that the solutions to (1) are non-negative. Theorem 3 ([13]). Uniformly bounded solutions for the model (1) begin in Ω = { (S,C, I,R) ∈ Ω+ : 0 ≤ S + C + I +R ≤ Λ µ } . Proof. Let N(t) = S(t) + C(t) + I(t) +R(t) be the total population at time t. Thus, we have DαN(t) = Λ− µN(t)− σI(t) ≤ Λ− µN(t). Thus, DαN(t) + µN(t) ≤ Λ. It follows from [18] that if Mα is the Mittag-Leffler function, we have 0 ≤ N(t) ≤ N(0)Mα(−µtα) + tαMα,α+1(−µtα)). Boukhouima et al. [15] yields 0 ≤ N(t) ≤ Λ µ , t −→ ∞. For this reason, the solutions of (1) are uniformly limited in the region Ω, beginning with Ω+. 3. Stability Analysis of the Model 3.1. Equilibria For I = C = 0, we have 0 = Λ− µS(t) + ηR(t), 0 = −(µ+ η)R(t). Therefore, S0 = Λ µ , R0 = 0. The disease-free equilibrium E1 = ( Λ µ , 0, 0, 0 ) is easy to obtain. The pneumonia endemic equilibrium is represented by E2 = (S∗, C∗, I∗, R∗) for I∗ > 0, and can be derived by solving the system of equations (1), that is, DαS(t) = 0, DαC(t) = 0, DαI(t) = 0, DαR(t) = 0. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2769 By taking α1 = θπ + (1− θ)Z1 Z2 , we have E∗ = (S∗, C∗, I∗, R∗), where S∗ = NZ1 α1 + ω , C∗ = (µ+ η)(Λ(α1 + ω)− µNZ1) (α1 + ω)(δZ1(µ+ η)− η(β + τα1)) , I∗ = α1(µ+ η)(Λ(α1 + ω)− µNZ1) (α1 + ω)(δZ1(µ+ η)− η(β + τα1)) , R∗ = (β + τα1)(Λ(α1 + ω)− µNZ1) (α1 + ω)(δZ1(µ+ η)− η(β + τα1)) . 3.2. The Basic Reproductive Number R0 Next-generation matrix theory allows us to deduce that we obtain Dαy = H(y)− h̄(y), if we set x = (S,R)T and y = (C, I)T , where H(y) = [ δθI(t)S(t) N + δωθC(t)S(t) N δ(I(t)+ωC(t)) N (1− θ)S(t) ] , h̄(y) = [ Z1C(t) −πC(t) + Z2I(t) ] . At E1 = ( Λ µ , 0, 0, 0 ) , the Jacobian matrix of H(y) and h̄(y) w.r.t. C and I may be obtained as follows: H(y) = [ H1 H2 ] , h̄(y) = [ h̄1 h̄2 ] , where H1 = δθI(t)S(t) N + δωθC(t)S(t) N , h̄1 = Z1C(t), H2 = δ(I(t) + ωC(t)) N (1− θ)S(t), h̄2 = −πC(t) + Z2I(t). Hence, at E0, the Jacobian matrix of H(y) and h̄(y) is derived by F and Γ, respectively, F = [ ∂H1 ∂C ∂H1 ∂I ∂H2 ∂C ∂H2 ∂I ] , Γ = [ ∂h̄1 ∂C ∂h̄1 ∂I ∂h̄2 ∂C ∂h̄2 ∂I ] . Then, ∂H1 ∂C = δωθ N S(t), ∂H1 ∂I = δθ N S(t), ∂H2 ∂C = δω(1− θ) N S(t), ∂H2 ∂I = δ(1− θ) N S(t). A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2770 Similarly, we may obtain the Γ entries as: ∂h̄1 ∂C = Z1, ∂h̄1 ∂I = 0, ∂h̄2 ∂C = −π, ∂h̄2 ∂I = Z2. The entries in A are derived as follows: A = [ δωθS0 N δθS0 N δω(1−θ)S0 N δ(1−θ)S0 N ] , B = [ Z1 0 −π Z2 ] . Thus, B−1 = [ 1 Z1 π Z1Z2 0 1 Z2 ] . Then, A.B−1 = [ δωθS0 N δθS0 N δω(1−θ)S0 N δ(1−θ)S0 N ] . [ 1 Z1 π Z1Z2 0 1 Z2 ] . A.B−1 = [ ΩθδS0 NZ1 (πω+Z1)θδS0 NZ1Z2 δω(1−θ)S0 NZ1 (πω+Z1)δ(1−θ)S0 NZ1Z2 ] . The eigenvalues of A.B−1 can be obtained as:∣∣∣∣∣λ− ΩθδS0 NZ1 (πω+Z1)θδS0 NZ1Z2 δω(1−θ)S0 NZ1 λ− (πω+Z1)δ(1−θ)S0 NZ1Z2 ∣∣∣∣∣ = 0. As a result, the eigenvalues are: λ1 = 0, λ2 = (ωθZ2 + (1− θ)(πω + Z1))δS0 NZ1Z2 . Thus, R0 = (ωθZ2 + (1− θ)(πω + Z1))δS0 NZ1Z2 . 3.3. Local Stability The Jacobian matrix of the system (1) along with its elements is presented as: J(En) =  J11 J12 J13 J14 J21 J22 J23 J24 J31 J32 J33 J34 J41 J42 J43 J44  , A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2771 where J11 = −δI(t) N − δωC(t) N − µ, J12 = −δωS(t) N , J13 = −δS(t) N , J14 = η, J21 = δθI(t) N + δωθC(t) N , J22 = δωθS(t) N − Z1, J23 = δθS(t) N , J24 = 0, J31 = δ(I(t) + ωC(t)) N (1− θ), J32 = δΩ N (1− θ)S(t) + π, J33 = δ N (1− θ)S(t)− Z2, J34 = 0, J41 = 0, J42 = β, J43 = τ, J44 = −(µ+ η). Lemma 2. The trivial disease equilibrium E0 = (0, 0, 0, 0) is locally asymptotically stable in Ω. Proof. At E0 = (0, 0, 0, 0), the model (1) has the following Jacobian matrix J(E0): J(E0) =  −µ 0 0 η 0 −Z1 0 0 0 π −Z2 0 0 β τ −(µ+ η)  . Therefore, we can get the eigenvalues as: J(E0) = ∣∣∣∣∣∣∣∣ λ+ µ 0 0 η 0 λ+ Z1 0 0 0 π λ+ Z2 0 0 β τ λ+ (µ+ η) ∣∣∣∣∣∣∣∣ = 0. Then, (λ+ µ)(λ+ µ+ η) ∣∣∣ λ+Z1 0 π λ+Z2 ∣∣∣ = 0. Therefore, (λ+ µ)(λ+ µ+ η)(λ+ Z1)(λ+ Z2) = 0. (2) From Equation (2), we have λ1 = −(λ+ µ), λ2 = −(λ+ µ+ η), λ3 = −Z1, λ4 = −Z2. Consequently, λ1, λ2, λ3, λ4 are all strictly negative roots of Eq. (2). Thus, as determined by the Routh-Hurwitz criterion, it follows that E0 = (0, 0, 0, 0) is locally asymptotically stable. Lemma 3. The disease-free equilibrium E1 = ( Λ µ , 0, 0, 0 ) is locally asymptotically stable in Ω if R0 < 1 and is unstable if R0 > 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2772 Proof. At E1 = ( Λ µ , 0, 0, 0 ) , the model (1) has the following Jacobian matrix J(E1): J(E1) = −µ − δωS0 N − δS0 N η 0 δωθS0 N −Z1 δθS0 N 0 0 δω(1−θ)S0 N +π δ(1−θ)S0 N −Z2 0 0 β τ −(µ+η)  . Therefore, we can get the eigenvalues as: J(E1) = ∣∣∣∣∣∣∣ λ+µ − δωS0 N − δS0 N η 0 λ− δωθS0 N +Z1 δθS0 N 0 0 δω(1−θ)S0 N +π λ− δ(1−θ)S0 N +Z2 0 0 β τ λ+(µ+η) ∣∣∣∣∣∣∣ = 0. Then, (λ+ µ)(λ+ µ+ η) ∣∣∣∣ λ− δωθS0 N +Z1 δθS0 N δω(1−θ)S0 N +π λ− δ(1−θ)S0 N +Z2 ∣∣∣∣ = 0. Therefore, (λ+ µ)(λ+ µ+ η) = 0, (3) or ∣∣∣∣ λ− δωθS0 N +Z1 δθS0 N δω(1−θ)S0 N +π λ− δ(1−θ)S0 N +Z2 ∣∣∣∣ = 0. From Eq. (3), we have λ1 = −(λ+ µ), λ2 = −(λ+ µ+ η), or λ2 + c1λ+ c2 = 0, where c1 = −δωθS0 N + Z1 − δ(1− θ)S0 N + Z2, c2 = −πδθS0 N − δωθZ2S0 N − δ(1− θ)Z1S0 N + Z1Z2, c2 = − [θπ + ωθZ2 + (1− θ)Z1]δS0 N + Z1Z2, where Z1 = µ+ β + π, Z2 = τ + µ+ σ, S0 = Λ µ . As a result, if c1 > 0 and c2 > 0, the Routh-Hurwitz criteria states that E1 is locally asymptotically stable when R0 < 1. When R0 < 1 holds, we must demonstrate that c1 > 0 and c2 > 0. R0 = (ωθZ2 + (1− θ)(πω + Z1))δS0 NZ1Z2 . A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2773 If R0 < 1, we have (ωθZ2 + Z1(1− θ))δS0 N < Z1Z2. If Γ = min{Z2, Z1}, we have (ωθ + (1− θ))δS0 N < Z1Z2 Γ . Thus, (ωθ + (1− θ))δS0 N < Z1 + Z2. Then, c1 = Z1 + Z2 − (ωθ + (1− θ))δS0 N > 0. If R0 < 1, we have (ωθZ2 + (πω + Z1)(1− θ))δS0 NZ1Z2 < 1, which implies (ωθZ2 + (πω + Z1)(1− θ))δS0 < NZ1Z2. That is, NZ1Z2 − [ωθZ2 + (πω + Z1)(1− θ)]δS0 > 0. That is, Z1Z2 − [πω(1− θ) + ωθZ2 + Z1(1− θ)]δS0 N > 0, where Z1 = µ+ β + π, Z2 = τ + µ+ σ, S0 = Λ µ . Thus, c2 > 0. For this reason, when R0 < 1, E1 is locally asymptotically stable and unstable when R0 > 1. Lemma 4. If R0 > 1, the pneumonia endemic equilibrium E2 = (S∗, C∗, I∗, R∗) is locally asymptotically stable in Ω. Proof. At E2 = (S∗, C∗, I∗, R∗), the model (1) has the following Jacobian matrix J(E2): J(E2) =  −f1 − µ −f2 −f3 η f1θ f2θ − Z1 f3θ 0 f1(1− θ) f2(1− θ) + π f3(1− θ)− Z2 0 0 β τ −(µ+ η)  , A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2774 where f1 = δI∗ N + δωC∗ N , f2 = δωS∗ N , f3 = δS∗ N , and its characteristic equation is given by: |J(E2)− λI| = 0. That is,  −f1 − µ− λ −f2 −f3 η f1θ f2θ − Z1 − λ f3θ 0 f1(1− θ) f2(1− θ) + π f3(1− θ)− Z2 − λ 0 0 β τ −(µ+ η)− λ  = η  f1θ f2θ − Z1 − λ f3θ f1(1− θ) f2(1− θ) + π f3(1− θ)− Z2 − λ 0 β τ  − (µ+ η + λ)  −f1 − µ− λ −f2 −f3 f1θ f2θ − Z1 − λ f3θ f1(1− θ) f2(1− θ) + π f3(1− θ)− Z2 − λ  = 0. Thus, the characteristic values are given by: λ4 + ξ1λ 3 + ξ2λ 2 + ξ3λ+ ξ4 = 0, where ξ1 = Z1 + Z2 − f2θ − f3(1− θ) + 2µ+ η + f1, ξ2 = Z1Z2 − f2Z2θ − f3Z1(1− θ)− f3πθ + (f1 + µ)(µ+ η) + f1f2θ + f1f3(1− θ), ξ3 = (2µ+ η + f1)(Z1Z2 − f2Z2θ − f3Z1(1− θ)− f3πθ) + (µ+ f1)(µ+ η)(Z1 + Z2 − f2θ) + (µ+ η)f1f2θ + f1f2θZ2 + (µ+ η)(f1f3(1− θ) + f1f3πθ + f1f3Z1(1− θ)), ξ4 = (µ+ f1)(µ+ η)(Z1Z2 − f2Z2θ − f3Z1(1− θ)− f3πθ) + f1f2θ(µ+ η)Z2 + (µ+ η)(f1f3πθ + f1f3Z1(1− θ)). According to the Routh-Hurwitz criterion for the fourth-degree polynomial, the provided constraint has been verified if R0 > 1. 3.4. Global Stability Theorem 4. The disease-free equilibrium E1 = ( Λ µ , 0, 0, 0 ) is globally asymptotically stable in Ω. Proof. Take the Lyapunov function with positive definiteness L2(S,C, I,R) = ( S − S0 − S0 ln S S0 ) . A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2775 From [58], we have DαL2(S,C, I,R) ≤ (S − S0 S ) DαS = (S − S0 S )( Λ− δI(t)S(t) N − δωC(t)S(t) N − µS(t) + ηR(t) ) . At E1 = ( Λ µ , 0, 0, 0 ) , we have DαL2(S, V, C, I, R) ≤ (S − S0 S ) DαS = (S − S0 S )( Λ− δ(I(t) + ωC(t)) N S(t)− µS(t) + ηR(t) ) = ( S − S0 )(Λ S − δ(I(t) + ωC(t)) N − µ+ ηR(t) S ) = ( S − S0 )(Λ S + ηR(t) S − Λ S0 − ηR(t) S0 ) = ( S − S0 )( − Λ SS0 ( S − S0 ) − ηR(t) SS0 ( S − S0 )) = − Λ SS0 ( S − S0 )2 − ηR SS0 ( S − S0 )2 . Thus, for all (S,C, I,R) ∈ Ω, DαL2(S,C, I,R) < 0. Therefore, it follows from [30] that E1 is globally asymptotically stable in Ω. Theorem 5. The pneumonia endemic equilibrium E2 = (S∗, C∗, I∗, R∗) is globally asymp- totically stable in Ω. Proof. Take the Lyapunov function with positive definiteness L2(S,C, I,R) = ( S − S∗ − S∗ ln S S∗ ) + ( C − C∗ − C∗ ln C C∗ ) + ( I − I∗ − I∗ ln I I∗ ) + ( R−R∗ −R∗ ln R R∗ ) . From [58], we have DαL2(S, V, C, I, R) ≤ (S − S∗ S ) DαS + (C − C∗ C ) DαC + (I − I∗ I ) DαI + (R−R∗ R ) DαR = (S − S∗ S )( Λ− δ(I + ωC) N S − µS + ηR ) + (C − C∗ C )(δ(I(t) + ωC) N θS − Z1C ) + (I − I∗ I )(δ(I + ωC) N (1− θ)S + πC − Z2I ) + (R−R∗ R )( βC + τI − (µ+ η)R ) A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2776 = ( S − S∗ )(Λ S − δ(I + ωC) N − µ+ ηR S ) + ( C − C∗ )(δ(IS) NC θ + δωS N θ − Z1 ) + ( I − I∗ )(δS N (1− θ) + δΩSC NI (1− θ) + πC I − Z2 ) + ( R−R∗ )(βC R + τI(t) R − (µ+ η) ) = ( S − S∗ )( Λ NS + ηR NS − Λ NS∗ − ηR NS∗ ) + ( C − C∗ )(δ(IS) NC θ − δ(IS) NC∗ θ ) + ( I − I∗ )(δΩSC NI (1− θ) + πC I − δΩSC NI∗ (1− θ)− πC I∗ ) + ( R−R∗ )(βC R + τI R − βC R∗ − τI R∗ ) = ( S − S∗ )( − Λ NSS∗ ( S − S∗ ) − ηR NSS∗ ( S − S∗ )) + ( C − C∗ )( − δ(ISθ) NCC∗ ( C − C∗ )) + ( I − I∗ )( − δω(1− θ)SC NII∗ ( I − I∗ ) − πC II∗ ( I − I∗ )) + ( R−R∗ )( − βC NRR∗ ( R−R∗ ) − τI NRR∗ ( R−R∗ )) = − Λ NSS∗ ( S − S∗ )2 − ηR NSS∗ ( S − S∗ )2 − δ(ISθ) NCC∗ ( C − C∗ )2 − δω(1− θ)SC NII∗ ( I − I∗ )2 − πC II∗ ( I − I∗ )2 − βC RR∗ ( R−R∗ )2 − τI RR∗ ( R−R∗ )2 . Thus, for all (S,C, I,R) ∈ Ω, DαL2(S,C, I,R) < 0. Therefore, it follows from [30] that E2 is globally asymptotically stable in Ω. 4. Numerical Computational Methods of Solution In this section, we describe an implementation of the numerical methods we will apply for solving the proposed model (1). The fractional integral Jα mf(t) of the function f : R+ −→ R is defined as: Jα mf(t) = 1 Γ(α) ∫ t m (t− θ)α−1f(θ)dθ, t ≥ m, A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2777 where α ∈ R+ is a non-integer order and Γ(z) = ∫ ∞ 0 e−ttz−1dt is the Euler Gamma function. It follows from [48] that the Caputo fractional derivative Dαf(t) of order α > 0, n− 1 < α < n, n ∈ N is defined as: c mDαf(t) = 1 Γ(n− α) ∫ t m f (n)(θ) (t− θ)α+1−n dθ, t ≥ m. 4.1. The Adams-Bashforth-Moulton Method (ABMM) We consider the problem, Dα x(t) = f(t, x(t)), t ∈ [0, T ], 0 < α ≤ 1, x(j)(0) = x (j) 0 , j = 0, 1, 2, ...., n− 1. (4) Let h = T k be the step size, where k is a positive integer and T > 0, and xj be the approximate solution of x(tj) at t = tj , where tj = jh, j = 0, 1, ..., k. Take into consideration this integral part, In+1 = ∫ tn+1 0 (tn+1 − θ)α−1g(θ) dθ, n = 0, 1, 2, ..., k − 1. This can be approximately estimated using the following method In+1 ≈ ∫ tn+1 0 (tn+1 − θ)α−1g̃n+1(θ) dθ, n = 0, 1, 2, ..., k − 1, where g̃n+1(θ) is the approximation of g(θ) on the interval [0, tn+1]. Diethelm et al. [20] carried out the preliminary study on the fractional Adams tech- nique for (1), which may be formulated in the following way: xpn+1 = k−1∑ j=0 tjn+1 j! x (j) 0 + hα Γ(α+ 2) n∑ j=0 ℓj,n+1f(tj , xj). xn+1 = k−1∑ j=0 tjn+1 j! y (j) 0 + hα Γ(α+ 2)  n∑ j=0 ℓj,n+1f(tj , xj) + ℓn+1,n+1f(tn+1, x p n+1)  . (5) Here ℓj,n+1 are given by ℓj,n+1 =  nα+1 − (n− α)(n+ 1)α, if j = 0, (n− j + 2)α+1 + (n− j)α+1 − 2(n− j + 1)α+1, if 1 ≤ j ≤ n, 1, if j = n+ 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2778 By taking Z1 = µ+ β + π, Z2 = τ + µ+ σ, (6) the model (1) may be reformulated in the sense of the ABMM, according to Eq. (5), as follows: Sn+1 = S0 + hα Γ(α+ 2) ( Λ− δ(Ipn+1 + ωCp n+1) N Sp n+1 − µSp n+1 + ηRp n+1 ) + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 ( Λ− δ(I(tj) + ωC(tj)) N S(tj)− µS(tj) + ηR(tj) ) , Cn+1 = C0 + hα Γ(α+ 2) (δ(Ipn+1 + ωCp n+1) N θSp n+1 − Z1C p n+1 ) + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 (δ(I(tj) + ωC(tj)) N θS(tj)− Z1C(tj) ) , In+1 = I0 + hα Γ(α+ 2) (δ(Ipn+1 + ωCp n+1) N (1− θ)Sp n+1 + πCp n+1 − Z2I p n+1 ) + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 (δ(I(tj) + ωC(tj)) N (1− θ)S(tj) + πC(tj)− Z2I(tj) ) , Rn+1 = R0 + hα Γ(α+ 2) ( βCp n+1 + τIpn+1 − (µ+ η)Rp n+1 ) + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 ( βC(tj) + τI(tj)− (µ+ η)R(tj) ) , where Sp n+1, C p n+1, I p n+1, and Rp n+1 are given below: Sp n+1 = S0 + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 ( Λ− δ(I(tj) + ωC(tj)) N S(tj)− µS(tj) + ηR(tj) ) , Cp n+1 = C0 + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 (δ(I(tj) + ωC(tj)) N θS(tj)− Z1C(tj) ) , Ipn+1 = I0 + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 (δ(I(tj) + ωC(tj)) N (1− θ)S(tj) + πC(tj)− Z2I(tj) ) , Rp n+1 = R0 + hα Γ(α+ 2) n∑ j=0 Ψj,n+1 ( βC(tj) + τI(tj)− (µ+ η)R(tj) ) . A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2779 4.2. The Generalized Euler Method (GEM) The extended Euler’s method for numerically solving the initial value problem (4) with the Caputo derivatives was derived by Odibat and Momani [43]. First, let’s look at the starting value problem. If x, Dαx, and D2αx are continuous on [0, b] and a1, . . . , an are positive integer constants, then it follows from (4) that for each value of t, there is a value a1, which implies x(t) = x(t0) + (Dα x(t))(t0) tα Γ(α+ 1) + (D2α x(t))(a1) t2α Γ(2α+ 1) . (7) Substituting (Dα x(t))(t0) = f(t0, x(t0)) and h = t1 into Eq. (7) yields an equation for x(t1): x(t1) = x(t0) + f(t0, x(t0)) hα Γ(α+ 1) + (D2α x(t))(a1) h2α Γ(2α+ 1) . If the step size h is sufficiently small, the second-order term (involving h2α) may be neglected, allowing us to get x(t1) = x(t0) + hα Γ(α+ 1) f(t0, x(t0)). Repeating this procedure yields a sequence of points that estimates the solution x(t). This iterative general formula for the GEM when tj+1 = tj + h is x(tj+1) = x(tj) + hα Γ(α+ 1) f(tj , x(tj)), (8) for j = 0, 1, ..., k− 1. The model (1) may be reformulated in the sense of GEM, according to Eq. (8), as follows: S(tj+1) = S(tj) + hα Γ(α+ 1) ( Λ− δ(I(tj) + ωC(tj)) N S(tj)− µS(tj) + ηR(tj) ) , C(tj+1) = C(tj) + hα Γ(α+ 1) (δ(I(tj) + ωC(tj)) N θS(tj)− Z1C(tj) ) , I(tj+1) = I(tj) + hα Γ(α+ 1) (δ(I(tj) + ωC(tj)) N (1− θ)S(tj) + πC(tj)− Z2I(tj) ) , R(tj+1) = R(tj) + hα Γ(α+ 1) ( βC(tj) + τI(tj)− (µ+ η)R(tj) ) . Through this, we have Sn+1 = Sn + hα Γ(α+ 1) hα Γ(α+ 1) ( Λ− δ(In + ωCn) N Sn − µSn + ηRn ) , Cn+1 = Cn + hα Γ(α+ 1) (δ(In + ωCn) N θSn − Z1Cn ) , In+1 = In + hα Γ(α+ 1) (δ(In + ωCn) N (1− θ)Sn + πCn − Z2In ) , Rn+1 = Rn + hα Γ(α+ 1) ( βCn + τIn − (µ+ η)Rn ) , A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2780 where 0 < hα Γ(α+1) < 1 is the step size and subject to the initial conditions, S(0) = S0, C(0) = C0, I(0) = I0, R(0) = R0. We have Sn+1 = Sn = S, Cn+1 = Cn = C, In+1 = In = I, and Rn+1 = Rn = R at a fixed point. 4.3. The Generalized Runge-Kutta Method (GRKM) The numerical scheme of Equation (4) of the GRKM takes the following form: xn+1 = xn + 1 6 ( K1 + 2K2 + 2K3 +K4 ) , (9) K1 = ℏf(tn, xn), K2 = ℏf(tn + 1 2 ℏ, xn + 1 2 K1), K3 = ℏf(tn + 1 2 ℏ, xn + 1 2 K2), K4 = ℏf(tn + ℏ, xn +K3), where h̄ = hα Γ(α+ 1) . In the sense of (9), the numerical scheme of (1) of the GRKM takes the following form: Sn+1 = Sn + 1 6 ( K1 + 2K2 + 2K3 +K4 ) , Cn+1 = Cn + 1 6 ( Q1 + 2Q2 + 2Q3 +Q4 ) , In+1 = In + 1 6 ( O1 + 2O2 + 2O3 +O4 ) , Rn+1 = Rn + 1 6 ( P1 + 2P2 + 2P3 + P4 ) , where K1, Q1, O1, P1,K2, Q2, O2, P2,K3, Q3, O3, P3,K4, Q4, O4, and P4 are given below: K1 = h̄f1(tn, Sn, Cn, In, Rn) = h̄ ( Λ− δ(In + ωCn) N Sn − µSn + ηRn ) , Q1 = h̄f2(tn, Sn, Cn, In, Rn) = h̄ (δ(In + ωCn) N θSn − Z1Cn ) , O1 = h̄f3(tn, Sn, Cn, In, Rn) = h̄ (δ(In + ωCn) N (1− θ)Sn + πCn − Z2In ) , P1 = h̄f4(tn, Sn, Cn, In, Rn) = h̄ ( βCn + τIn − (µ+ η)Rn ) , K2 = h̄f1 ( tn + 1 2 h̄, Sn + 1 2 K1, Cn + 1 2 Q1, In + 1 2 O1, Rn + 1 2 P1 ) , Q2 = h̄f2 ( tn + 1 2 h̄, Sn + 1 2 K1, Cn + 1 2 Q1, In + 1 2 O1, Rn + 1 2 P1 ) , O2 = h̄f3 ( tn + 1 2 h̄, Sn + 1 2 K1, Cn + 1 2 Q1, In + 1 2 O1, Rn + 1 2 P1 ) , P2 = h̄f4 ( tn + 1 2 h̄, Sn + 1 2 K1, Cn + 1 2 Q1, In + 1 2 O1, Rn + 1 2 P1 ) , A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2781 K3 = h̄f1 ( tn + 1 2 h̄, Sn + 1 2 K2, Cn + 1 2 Q2, In + 1 2 O2, Rn + 1 2 P2 ) , Q3 = h̄f2 ( tn + 1 2 h̄, Sn + 1 2 K2, Cn + 1 2 Q2, In + 1 2 O2, Rn + 1 2 P2 ) , O3 = h̄f3 ( tn + 1 2 h̄, Sn + 1 2 K2, Cn + 1 2 Q2, In + 1 2 O2, Rn + 1 2 P2 ) , P3 = h̄f4 ( tn + 1 2 h̄, Sn + 1 2 K2, Cn + 1 2 Q2, In + 1 2 O2, Rn + 1 2 P2 ) , K4 = h̄f1(tn + h̄, Sn +K3, Cn +Q3, In +O3, Rn + P3), Q4 = h̄f2(tn + h̄, Sn +K3, Cn +Q3, In +O3, Rn + P3), O4 = h̄f3(tn + h̄, Sn +K3, Cn +Q3, In +O3, Rn + P3), P4 = h̄f4(tn + h̄, Sn +K3, Cn +Q3, In +O3, Rn + P3). 4.4. The Multistep Generalized Differential TransformMethod (MSGDTM) To illustrate the MSGDTM following ([42], [24], [22], [25]), we consider the following initial value problem for systems of fractional differential equations: Dα 1 x1(t) = f1(t, x1, x2, ..., xn), Dα 2 x1(t) = f2(t, x1, x2, ..., xn), ... Dα n x1(t) = fn(t, x1, x2, ..., xn), (10) subject to the initial conditions: xi(t0) = ci, i = 1, 2, ..., n, (11) where Dα is the Caputo fractional derivative of order αi, i = 1, 2, ..., n. The initial value problem solution (10), (11) is looked in the range [t0, T ]. In real-world applications of the extended differential transform technique, the finite series defines the kth order approxi- mation solution of (4) xi(t) = K∑ i=0 Li(k)(t− t0) kαi , t ∈ [t0, T ], where Li(k) satisfied the (recurrence relation): Γ((k + 1)αi + 1) Γ(kαi + 1) Li(k + 1) = Fi(k,L1,L2, ...,Ln), Li(0) = ci, and Fi(k,L1,L2, ...,Ln) is the differential transform of function fit, fi(t, x1, x2, ..., xn) for i = 1, 2, ..., n. The basic steps of the GDTM can be found in ([41], [37], [40], [23]). A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2782 Assume that [t0, T ] is partitioned into M subintervals [tm−1, tm], m = 1, 2, ...,M , with steps h = T − t0 M using the nodes tm = t0 + mh. The MSGDTM principles are listed below. First, the GDTM solves (10), (11) throughout the interval [t0, t1]. We get the approxi- mate solution by beginning with xi,1(t), t ∈ [t0, t1], for i = 1, 2, ..., n. For m ≥ 2, applying the GDTM to (10), (11) across the interval [tm−1, tm] will employ the IC xi,m(tm−1) = xi,m−1(tm−1). Repeating the technique yields approximations for xi,m(t), m = 1, 2, ...,M , for i = 1, 2, ..., n. The MSGDTM concludes: xi(t) =  xi,1(t), for t ∈ [t0, t1], xi,2(t), for t ∈ [t1, t2], ... xi,M (t), for t ∈ [tM−1, tM ]. The MSGDTM is efficient for all values of h. The new algorithm’s key benefit is that the solution converges over vast time periods, as we shall show in the following section. Regarding (1), the MSGDTM algorithm produces S(k + 1) = Γ(q1k + 1) Γ(q1(k + 1) + 1) ( Λ− δ(I(k) + ωC(k)) N S(k)− µS(k) + ηR(k) ) , C(k + 1) = Γ(q2k + 1) Γ(q2(k + 1) + 1) (δ(I(k) + ωC(k)) N θS(k)− Z1C(k) ) , I(k + 1) = Γ(q3k + 1) Γ(q3(k + 1) + 1) (δ(I(k) + ωC(k)) N (1− θ)S(k) + πC(k)− Z2I(k) ) , R(k + 1) = Γ(q4k + 1) Γ(q4(k + 1) + 1) ( βC(k) + τI(k)− (µ+ η)R(k) ) , where S(k), C(k), I(k), and R(k) are the differential transformations of S(k), C(k), I(k), and R(k), respectively. The differential transform of the starting conditions is provided by the formulas S(0) = c1, C(0) = c2, I(0) = c3, and R(0) = c4, respectively. Through the differential inverse transform, we may find the differential transform series solution for the problem (1): S(t) = K∑ n=0 S(n)tq1n, C(t) = K∑ n=0 C(n)tq2n, I(t) = K∑ n=0 I(n)tq3n, R(t) = K∑ n=0 R(n)tq4n. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2783 For the MSGDTM, the proposed series solution of (1) is: S(t) =  ∑K n=0 S1(n)tq1n, for t ∈ [0, t1]∑K n=0 S2(n)(t− t1) q1n, for t ∈ [t1, t2] ...∑K n=0 SM (n)(t− tM−1) q1n, for t ∈ [tM−1, tM ] (12) C(t) =  ∑K n=0C1(n)t q2n, for t ∈ [0, t1]∑K n=0C2(n)(t− t1) q2n, for t ∈ [t1, t2] ...∑K n=0CM (n)(t− tM−1) q2n, for t ∈ [tM−1, tM ] (13) I(t) =  ∑K n=0 I1(n)tq3n, for t ∈ [0, t1]∑K n=0 I2(n)(t− t1) q3n, for t ∈ [t1, t2] ...∑K n=0 IM (n)(t− tM−1) q3n, for t ∈ [tM−1, tM ] (14) R(t) =  ∑K n=0R1(n)t q4n, for t ∈ [0, t1]∑K n=0R2(n)(t− t1) q4n, for t ∈ [t1, t2] ...∑K n=0RM (n)(t− tM−1) q4n, for t ∈ [tM−1, tM ], (15) where Si(n), Ci(n), Ii(n), and Ri(n) for i = 1, 2, ...,M satisfy the following (recurrence relations): Si(k + 1) = Γ(q1k + 1) Γ(q1(k + 1) + 1) ( Λ− δ(Ii(k) + ωCi(k)) N Si(k)− µSi(k) + ηRi(k) ) , Ci(k + 1) = Γ(q2k + 1) Γ(q2(k + 1) + 1) (δ(Ii(k) + ωCi(k)) N θSi(k)− Z1Ci(k) ) , Ii(k + 1) = Γ(q3k + 1) Γ(q3(k + 1) + 1) (δ(Ii(k) + ωCi(k)) N (1− θ)Si(k) + πCi(k)− Z2Ii(k) ) , Ri(k + 1) = Γ(q4k + 1) Γ(q4(k + 1) + 1) ( βCi(k) + τIi(k)− (µ+ η)Ri(k) ) , (16) with Si(0) = Si(ti−1) = Si−1(ti−1), Ci(0) = Ci(ti−1) = Ci−1(ti−1), Ii(0) = Ii(ti−1) = Ii−1(ti−1), and Ri(0) = Ri(ti−1) = Ri−1(ti−1). Finally, we start with S0(0) = c1, C0(0) = c2, I0(0) = c3, and R0(0) = c4, using the (recurrence relations) (16), then we can get the multistep solution (12)-(15). A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2784 5. Numerical Simulations and Results Using the mathematical software programs: MATLAB and Mathematica and the pa- rameter values listed in Table 1, which are taken from [1], we present numerical simulations and results in this section. Table 1: The Parameter Values Parameter Value / Source [1] Λ 0.5 δ 2 (DFE) 2.5 (EE) ω 0.1124 µ 0.5 η 0.00641 β 0.515 π 0.7096 θ 0.563 τ 0.641 σ 0.53 5.1. Numerical Simulations and Results of the ABMM Method In this subsection, we present numerical simulations and results for the model (1) solution applying the ABMM method. The solution to the proposed model (1) using the Adams-Bashforth-Moulton Method appears in Figures 1(a, b, c, d) and Figures 2(a, b, c, d, e). According to these figures, we observe that S(t) increases at α = 1 than in the case where α is a fractional variant. We also note that C(t), I(t), and R(t) decrease at α = 1 than in the case where α is a fractional variant. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2785 0 10 20 30 40 50 60 70 80 90 100 Time (days) 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 S (a) 0 10 20 30 40 50 60 70 80 90 100 Time (days) -0.02 0 0.02 0.04 0.06 0.08 0.1 C (b) 0 10 20 30 40 50 60 70 80 90 100 Time (days) -0.02 0 0.02 0.04 0.06 0.08 0.1 I (c) 0 10 20 30 40 50 60 70 80 90 100 Time (days) -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 R (d) Figure 1: The ABMM method for the behavior of S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2786 Time,α=0.65 0 10 20 30 40 50 60 S , C , I, R 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 S C I R (a) Time,α=0.75 0 10 20 30 40 50 60 S , C , I, R 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 S C I R (b) Time,α=0.85 0 10 20 30 40 50 60 S , C , I, R 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 S C I R (c) Time,α=0.95 0 10 20 30 40 50 60 S , C , I, R 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 S C I R (d) Time,α=1 0 10 20 30 40 50 60 S , C , I, R -0.2 0 0.2 0.4 0.6 0.8 1 S C I R (e) Figure 2: Combined graphical behaviors of the ABMM method for S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2787 5.2. Numerical Simulations and Results of the GEM Method In this subsection, we present numerical simulations and results for the model (1) solution applying the GEM method. The solution to the proposed model (1) using the Generalized Euler Method appears in Figures 3(a, b, c, d) and Figures 4(a, b, c, d, e). According to these figures, we observe that S(t) increases at α = 1 than in the case where α is a fractional variant. We also note that C(t) and I(t) decrease at α = 1 than in the case where α is a fractional variant. For R(t), the dynamical behavior increases at α = 1 in the first days and then gets remarkably decreased, as well as in the case where α is a fractional variant. 0 20 40 60 80 100 0.90 0.92 0.94 0.96 0.98 1.00 Time (days) S �=0.65 �=0.75 (a) 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 Time (days) C �=1 (b) 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 Time (days) I �=1 (c) 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Time (days) R �=0.65 �=0.75 �=0.85 �=0.95 �=1 (d) Figure 3: The GEM method for the behavior of S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2788 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 Time (days),�=0.65 S ,C ,I ,R R I C S (a) 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 Time (days),�=0.75 S ,C ,I ,R R I C S (b) 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 Time (days),�=0.85 S ,C ,I ,R R I C S (c) 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 Time (days),�=0.95 S ,C ,I ,R R I C S (d) 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 Time (days),�=1 S ,C ,I ,R R I C S (e) Figure 4: Combined graphical behaviors of the GEM method for S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2789 5.3. Numerical Simulations and Results of the GRKM Method In this subsection, we present numerical simulations and results for the model (1) solution applying the GRKM method. In Figures 5(a, b, c, d) and Figures 6(a, b, c, d, e), we present the solution to the proposed model (1) using the Generalized Runge-Kutta Method. According to these figures, we observe that S(t) increases more in the case where α is a fractional variant than in the case where α = 1. We also note that C(t), I(t), and R(t) decrease more in the case where α is a fractional variant than in the case where α = 1. (a) (b) (c) (d) Figure 5: The GRKM method for the behavior of S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2790 (a) (b) (c) (d) (e) Figure 6: Combined graphical behaviors of the GRKM method for S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2791 5.4. Numerical Simulations and Results of the MSGDTM Method In this subsection, we present numerical simulations and results for the model (1) solution applying the MSGDTM method. The solution to the proposed model (1) using the Multistep Generalized Differential TransformMethod appears in Figures 7(a, b, c, d) and Figures 8(a, b, c, d, e). According to these figures, we observe that S(t) remarkably increases in the case where α is a fractional variant than in the case where α = 1. We also note that C(t), I(t), and R(t) remarkably decrease in a similar way for both cases where α is a fractional variant and α = 1. (a) (b) (c) (d) Figure 7: The MSGDTM method for the behavior of S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2792 (a) (b) (c) (d) (e) Figure 8: Combined graphical behaviors of the MSGDTM method for S(t), C(t), I(t), and R(t) at various values of α = 0.65, 0.75, 0.85, 0.95, 1. A. Alalyani / Eur. J. Pure Appl. Math, 17 (4) (2024), 2763-2799 2793 6. Discussion of the Results We present the solution to the proposed model (1) using four different numerical meth- ods. From the numerical simulations and results in Section 5, we observe that the solution to the proposed model (1) gets stable faster using the ABMM and the MSGDTM methods than using the GEM and the GRKM methods. The solution to the proposed model (1) using the ABMM method gets stable faster in the case where α = 1 than in the case where α is a fractional variant. On the other hand, the solution to the proposed model (1) using the MSGDTM method gets stable faster in the case where α is a fractional variant than in the case where α = 1. From numerical simulations and results, we can say that the numerical methods we used are efficient computational methods. The ABMM, the GEM, the GRKM, and the MSGDTM exhibit greater flexibility and applicability across a diverse spectrum of problems, with a specific emphasis on initial value problems and systems of ordinary differential equations (ODEs). The improved shooting method, specialized for boundary value problems (BVPs), exhibits high accuracy, stability, and robustness when applied to coupled nonlinear higher-order systems. It is generally the preferred method for dealing with BVPs, particularly in complex nonlinear intricate. 7. Conclusion In this work, we investigated the solution and the dynamics of the pneumonia fractional- order mathematical model (1) using numerical techniques including the Adams-Bashforth- Moulton method, the generalized Euler method, the generalized Runge-Kutta method, and the multistep generalized differential transform method. For the sake of the analysis of the proposed model (1), we discussed positivity, boundedness, equilibria, local and global stability, and the basic reproductive number R0. Based on the numerical simulations, we observed that there is a convergence in results between the proposed fractional-order mathematical model (1) and its integer-order form. The aim of using the fractional-order mathematical model (1) is to investigate the memory of the model for a long time period; that is, at any time t, the fractional-order mathematical model investigates the behavior of the variables S(t), C(t), I(t), and R(t) during the period [0, 1] as a result of changing the values of α, while the integer-order form investigates the behavior of the variables only at α = 1. This was useful in investigating the stability of the model during the period [0, 1] as a result of changing the values of α and gaining a better understanding of the disease’s dynamics. This work will also be a valuable tool for different types of disease modeling. The figures presented in this paper illustrate the dynamic behavior of the system variables S(t), C(t), I(t), and R(t) under various fractional orders α using four different numerical methods: the Adams-Bashforth-Moulton Method (ABMM), the Generalized Euler Method (GEM), the Generalized Runge-Kutta Method (GRKM), and the Multistep Generalized Differential Transform Method (MSGDTM). The main contribution of these figures lies in conducting a comparative analysis of the four different methods used to approximate solutions for the proposed fractional-order mathematical model. Specifically, Figures 1 and 2, obtained using the ABMM method, REFERENCES 2794 and Figures 3 and 4, obtained employing the GEM method, illustrate that S(t) exhibits an increase, while C(t), I(t), and R(t) all decrease as the fractional order α approaches 1. This observation indicates that the behavior of integer-order dynamics diverges notably from that of fractional-order dynamics, and the depicted figures highlight the model’s sensitivity to fractional variations. Figures 5 and 6, obtained using the GRKM method, illustrate a reverse trend where S(t) increases more for fractional values of α compared to α = 1, while C(t), I(t), and R(t) decrease, emphasizing the impact of employing a fractional order in the model. Figures 7 and 8, obtained using the MSGDTM method, demonstrate a remarkable increase in S(t) for fractional values of α compared to α = 1, while C(t), I(t), and R(t) consistently decrease, reinforcing the consistent behavior observed across various methods. The collective results indicate that fractional-order models offer a more nuanced and accurate representation of system dynamics when compared to integer-order models. Additionally, the selection of a numerical method significantly impacts the accuracy and stability of the solutions. The visual comparisons presented in the figures play a pivotal role in comprehending distinctions and in choosing suitable numerical techniques for fractional- order differential equations. Acknowledgements The author expresses gratitude to the editor and referees for their time and insightful comments. References [1] K. Abodayeh, A. Raza, M. Rafiq, M. Shoaib Arif, M. Naveed, et al. Analysis of Pneu- monia Model via Efficient Computing Techniques. Computers, Materials & Continua, 70(3):6073–6088, 2022. [2] Khalid Ahmed, Haroon Adam, M.Y. Youssif, and Sayed Saber. Different strategies for diabetes by mathematical modeling: Applications of fractal-fractional derivatives in the sense of atangana-baleanu. Results in Physics, 52:106892, 2023. [3] Khalid Ahmed, Haroon Adam, M.Y. Youssif, and Sayed Saber. Different strategies for diabetes by mathematical modeling: Modified minimal model. Alex. Eng. J., 80:74–87, 2023. [4] Khalid I.A. Ahmed, Haroon D.S. Adam, Najat Almutairi, and Sayed Saber. Analyt- ical solutions for a class of variable-order fractional liu system under time-dependent variable coefficients. Results in Physics, 56:107311, 2024. [5] S. Al-Zahrani, Fathelrhman Guma, K. Al-Zahrani, and Sayed Saber. A Fractional Or- der SITR Model for Forecasting of Transmission of COVID-19: Sensitivity Statistical Analysis. Malaysian Journal of Mathematical Sciences, 16(3):517–536, 2022. REFERENCES 2795 [6] Ahmad Alalyani. On the Solution of a Nonlinear Fractional-Order Glucose-Insulin System Incorporating β-cells Compartment. Malaysian Journal of Mathematical Sci- ences, 17(1):1–12, 2023. [7] Ahmad Alalyani and Sayed Saber. Stability analysis and numerical simulations of the fractional COVID-19 pandemic model. Int. J. Nonlin. Sci. Num., 24(3):989– 1002, 2023. [8] N. Almutairi and S. Saber. On chaos control of nonlinear fractional newton-leipnik system via fractional caputo-fabrizio derivatives. Scientific Reports, 13:22726, 2023. [9] Najat Almutairi and Sayed Saber. Chaos control and numerical solution of time- varying fractional Newton-Leipnik system using fractional Atangana-Baleanu deriva- tives. AIMS Mathematics, 8(11):25863–25887, 2023. [10] Najat Almutairi and Sayed Saber. Application of a time-fractal fractional derivative with a power-law kernel to the burke-shaw system based on newton’s interpolation polynomials. MethodsX, 12:102510, 2024. [11] Najat Almutairi and Sayed Saber. Existence of chaos and the approximate solution of the Lorenz–Lü–Chen system with the Caputo fractional operator. AIP Advances, 14(1):015112, 2024. [12] Najat Almutairi, Sayed Saber, and Hijaz Ahmad. The fractal-fractional Atangana- Baleanu operator for pneumonia disease: stability, statistical and numerical analyses. AIMS Mathematics, 8(12):29382–29410, 2023. [13] M. H. Alshehri, F. Z. Duraihem, A. Alalyani, and S. Saber. A Caputo (discretiza- tion) fractional-order model of glucose-insulin interaction: Numerical solution and comparisons with experimental data. J. Taibah Univ. Sci., 15(1):26–36, 2021. [14] Mansoor H. Alshehri, Sayed Saber, and Faisal Z. Duraihem. Dynamical analysis of fractional-order of ivgtt glucose–insulin interaction. Int. J. Nonlin. Sci. Num., 24:1123–1140, 2023. [15] A. Boukhouima, K. Hattaf, and N. Yousfi. Dynamics of a fractional order hiv infection model with specific functional response and cure rate. Int. J. Differ. Equ., 2017:Article ID 8372140, 2017. [16] A.C.G. Cesar, L.F.C. Nascimento, K.C.C. Mantovani, and L.C.P. Vieira. Fine par- ticulate matter estimated by mathematical model and hospitalisations for pneumonia and asthma in children. Rev. Paul. Pediatr. (Engl. Ed.), 34(1):18–23, 2016. [17] Y.H. Cheng, S.H. You, Y.J. Lin, S.C. Chen, W.Y. Chen, W.C. Chou, and C.M. Liao. Mathematical modeling of post coinfection with influenza a virus and streptococ- cus pneumoniae, with implications for pneumonia and copd-risk assessment. Int. J. Chronic Obstr. Pulm. Dis., 12:1973–1988, 2017. REFERENCES 2796 [18] S.K. Choi, B. Kang, and N. Koo. Stability for caputo fractional differential systems. Abstr. Appl. Anal., 2014:Article ID 631419, 2014. [19] I.M. Diah and N. Aziz. Stochastic modelling for pneumonia incidence: a conceptual framework. AIP Conference Proceedings New York, 1:1–3, 2019. [20] K. Diethelm and A. D. Freed. The fracpece subroutine for the numerical solution of differential equations of fractional order. Forschung und wissenschaftliches Rechnen, 1999:57–71, 2002. [21] G.L. Drusano, W. Liu, S. Fikes, R. Cirz, N. Robbins, S. Kurhanewicz, and A. Louie. Interaction of drug-and granulocyte-mediated killing of pseudomonas aeruginosa in a murine pneumonia model. J. Infect. Dis., 210(8):1319–1324, 2014. [22] V. Erturk, G. Zaman, and S. Momani. A numeric-analytic method for approximating a giving up smoking model containing fractional derivatives. Computers & Mathe- matics with Applications, 64(10):3065–3074, 2012. [23] V. S. Erturk, S. Momani, and Z. Odibat. Application of generalized differential transform method to multi-order fractional differential equations. Communications in Nonlinear Science and Numerical Simulation, 13(8):1642–1654, 2008. [24] V. S. Erturk, Z. M. Odibat, and S. Momani. An approximate solution of a fractional order differential equation model of human T-cell lymphotropic virus I (HTLV-I) infection of CD4 T-cells. Computers & Mathematics with Applications, 62(3):996– 1002, 2011. [25] Asad Freihat and Shaher Momani. Application of multistep generalized differential transform method for the solutions of the fractional-order chua’s system. Discrete Dynamics in Nature and Society, 2012:Article ID 427393, 2012. [26] S.M. Jung, R. Kinoshita, R.N. Thompson, N.M. Linton, Y. Yang, A.R. Akhmet- zhanov, and H. Nishiura. Epidemiological identification of a novel pathogen in real- time: analysis of the atypical pneumonia outbreak in Wuhan, China, 2019-2020. J. Clin. Med., 9(3):1–18, 2020. [27] M. Kizito and J. Tumwiine. A mathematical model of treatment and vaccination inter- ventions of pneumococcal pneumonia infection dynamics. J. Appl. Math., 2018:Article ID 2539465, 2018. [28] K. Kosasih and U. Abeyratne. Exhaustive mathematical analysis of simple clinical measurements for childhood pneumonia diagnosis. World J. Pediatr., 13(5):446–456, 2017. [29] K. Kosasih, U.R. Abeyratne, V. Swarnkar, and R. Triasih. Wavelet augmented cough analysis for rapid childhood pneumonia diagnosis. IEEE Trans. Biomed. Eng., 62(4):1185–1194, 2015. REFERENCES 2797 [30] J.P. LaSalle. The Stability of Dynamical Systems. SIAM, Philadelphia, 1976. [31] H.L. Li, L. Zhang, C. Hu, et al. Dynamical analysis of a fractional-order predator- prey model incorporating a prey refuge. J. Appl. Math. Comput., 54(1-2):435–449, 2017. [32] C. Marchello, A.P. Dale, T.N. Thai, D.S. Han, and M.H. Ebell. Prevalence of atyp- ical pathogens in patients with cough and community-acquired pneumonia: a meta- analysis. Ann. Fam. Med., 14(6):552–566, 2016. [33] F.K. Mbabazi, J.Y.T. Mugisha, and M. Kimathi. Modeling the within-host coinfection of influenza a virus and pneumococcus. Appl. Math. Comput., 339:488–506, 2018. [34] Johns Hopkins Medicine. Johns hopkins medicine: Pneumonia. Report n.d., Johns Hopkins Medicine, n.d. [35] W.K. Ming, J. Huang, and C.J. Zhang. Breaking down of healthcare system: math- ematical modelling for controlling the novel coronavirus (2019-ncov) outbreak in wuhan, china. BioRxiv, 1:1–18, 2020. [36] E. Mochan, D. Swigon, G. Ermentrout, S. Luken, and G.A. Clermont. Mathematical model of intrahost pneumococcal pneumonia infection dynamics in murine strains. J. Theor. Biol., 353:44–54, 2014. [37] S. Momani and Z. Odibat. A novel method for nonlinear fractional partial differen- tial equations: combination of DTM and generalized Taylor’s formula. Journal of Computational and Applied Mathematics, 220(1-2):85–95, 2008. [38] M. Naveed, D. Baleanu, A. Raza, et al. Modeling the transmission dynamics of delayed pneumonia-like diseases with a sensitivity of parameters. Adv. Differ. Equ., 2021:468, 2021. [39] E.J. Ndelwa, M. Kgosimore, E.S. Massawe, and L. Namkinga. Mathematical mod- elling and analysis of treatment and screening of pneumonia. Math. Theory Model., 5(10):21–39, 2015. [40] Z. Odibat and S. Momani. A generalized differential transform method for linear par- tial differential equations of fractional order. Applied Mathematics Letters, 21(2):194– 199, 2008. [41] Z. Odibat, S. Momani, and V. S. Erturk. Generalized differential transform method: application to differential equations of fractional order. Applied Mathematics and Computation, 197(2):467–477, 2008. [42] Z. M. Odibat, C. Bertelle, M. A. Aziz-Alaoui, and G. H. E. Duchamp. A multi- step differential transform method and application to non-chaotic or chaotic systems. Computers & Mathematics with Applications, 59(4):1462–1472, 2010. REFERENCES 2798 [43] Z.M. Odibat and S. Momani. An algorithm for the numerical solution of differential equations of fractional order. J. Appl. Math. Inform., 26(1-2):15–27, 2008. [44] J.P. Olumuyiwa, A. Yusuf, K. Oshinubi, F.A. Oguntolu, J.O. Lawal, A.I. Abioye, and T.A. Ayoola. Fractional order of pneumococcal pneumonia infection model with caputo fabrizio operator. Results in Physics, 29:104–581, 2021. [45] K.I. Oluwatobi and L.M. Erinle-Ibrahim. Mathematical modeling of pneumonia dy- namics of children under the age of five. Res. Square, 1:1–16, 2021. [46] J. Ong’ala, P. Oleche, and J.Y.T. Mugisha. Mathematical model for pneumonia dynamics with carriers. Int. J. Math. Anal., 7(50):2457–2473, 2013. [47] D. Otoo, P. Opoku, S. Charles, and A.P. Kingsley. Deterministic epidemic model for (svcsycasyir) pneumonia dynamics, with vaccination and temporal immunity. Infect. Dis. Model., 5:42–60, 2020. [48] I. Podlubny. Fractional Differential Equations. Academic Press, New York, NY, USA, 1999. [49] M. Raj, M. Reddy, M. Mufeed, and S. Karthika. HMM based cough sound analysis for classifying asthma and pneumonia in paediatric population. Int. J. Pure Appl. Math., 118(18):609–616, 2018. [50] S. Saber and A. Alalyani. Stability analysis and numerical simulations of ivgtt glucose- insulin interaction models with two time delays. Math. Model. Anal., 27:383–407, 2022. [51] Sayed Saber. Control of Chaos in the Burke-Shaw system of fractal-fractional order in the sense of Caputo-Fabrizio. Journal of Applied Mathematics and Computational Mechanics, 23(1):83–96, 2024. [52] Sayed Saber, Azza M. Alghamdi, Ghada A. Ahmed, and Khulud M. Alshehri. Math- ematical modelling and optimal control of pneumonia disease in sheep and goats in Al-Baha region with cost-effective strategies. AIMS Mathematics, 7(7):12011–12049, 2022. [53] G.T. Tilahun. Modeling co-dynamics of pneumonia and meningitis diseases. Adv. Differ. Equ., 2019(1):1, 2019. [54] G.T. Tilahun. Optimal control analysis of pneumonia and meningitis coinfection. Comput. Math. Methods Med., 2019:1–15, 2019. [55] G.T. Tilahun, O.D. Makinde, and D. Malonza. Modelling and optimal control of pneumonia disease with cost-effective strategies. J. Biol. Dyn., 11(2):400–426, 2017. REFERENCES 2799 [56] G.T. Tilahun, O.D. Makinde, and D. Malonza. Co-dynamics of pneumonia and ty- phoid fever diseases with cost-effective optimal control analysis. Appl. Math. Comput., 316:438–459, 2018. [57] R. Ullah, R. Ellahi, S. M. Sait, and S. T. Mohyud-Din. On the fractional order model of HIV-1 infection of CD4+ T-cells under the influence of antiviral drug treatment. J. Taibah. Univ. Sci., 14(1):50–59, 2020. [58] C. Vargas-De-León. Volterra-type lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simul., 24:75–85, 2015. [59] N.M. Wafula, B.O. Kwach, and V.N. Marani. Mathematical modeling and opti- mal control for controlling pneumonia- HIV coinfection. Int. J. Innov. Res. Dev., 10(1):138–144, 2021. [60] O.C. Zephaniah, U.I.R. Nwaugonma, I.S. Chioma, and O. Adrew. A mathematical model and analysis of an SVEIR model for streptococcus pneumonia with saturated incidence force of infection. Math. Model. Appl., 5(1):16, 2020.