EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6176 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Hybrid Exponential Runge-Kutta Scheme for Stochastic SIQR Disease Modeling Muhammad Shoaib Arif1,∗, Kamaleldin Abodayeh1, Yasir Nawaz2 1 Department of Mathematics and Sciences, College of Humanities and Sciences, Prince Sultan University, Riyadh, 11586, Saudi Arabia 2 Department of Mathematics, Air University, PAF Complex E-9, Islamabad, 44000, Pakistan Abstract. Understanding and reducing the spread of epidemics depends much on the modelling and study of infectious disease dynamics. Among several compartmental models, the Susceptible- Infected-Quarantined-Recovered SIQR framework has become more popular as it can include the influence of quarantine, a significant intervention in many actual epidemics. Fundamental epidemic processes are naturally vulnerable to random variations brought on by environmental variability, population stochasticity, and other unknown elements. Hence, including stochastic influences in such models is crucial. Furthermore, spatial dispersion and diffusion effects are essential, partic- ularly in significant populations and varied settings, which call for stochastic partial differential equations (SPDEs). A two-stage mixture of exponential integrator and Runge-Kutta scheme is proposed for solving stochastic epidemic disease models. The scheme is more accurate than the existing Euler Maruyama method. The stability and consistency of the scheme in the mean square sense are provided. The scheme only discretizes time-dependent terms in given stochastic partial differential equations. Moreover, a stochastic SIQR diffusive model is presented with the effect of incidence rate. The deterministic and stochastic models are solved using the Euler-Maruyama method, the proposed scheme, and the nonstandard finite difference method. The comparison shows that the proposed scheme provides less error than the existing nonstandard finite difference method. The results indicate that the proposed scheme attains superior accuracy and diminished errors relative to current methods. This underscores its capability as an effective instrument for simulating complex stochastic epidemic models incorporating spatial effects and non-linear dynam- ics. 2020 Mathematics Subject Classifications: 65C30, 35R60, 92D30 Key Words and Phrases: Exponential integrator scheme, Runge-Kutta scheme, stability, mean square consistency, SIQR model, computational epidemiology ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6176 Email address: marif@psu.edu.sa (M. S. Arif), kamal@psu.edu.sa (K. Abodayeh), yasir maths@yahoo.com (Y. Nawaz) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 2 of 21 1. Introduction In recent decades, the increasing number of infectious diseases has created significant difficulties for public health systems. Mathematical modelling has greatly aided in un- derstanding the dynamics of disease transmission and assessing control measures. Among many compartmental models, the SIQR (Susceptible-Infected-Quarantined-Recovered) model has attracted significant interest because it can reflect the consequences of recovery pro- cesses and quarantine actions. Although insightful, conventional deterministic models can overlook the natural randomness in actual epidemic data. This constraint emphasizes the need for stochastic modelling, which considers random variations in disease transmission and other epidemiological factors and offers a more realistic picture of disease dynamics. Particularly in situations with small populations or unknown settings, the standard SIQR model’s incorporation of stochastic components improves the prediction and analysis of epidemic trends. Furthermore, the model’s relevance is increased by considering a generalized incidence rate instead of the straightforward bilinear form, capturing non- linear infection patterns seen in real epidemics. Accurately portraying situations where new infection rates are affected by variables, including behavioural changes, demographic diversity, or public health measures, requires this shift. The main objectives of this research are as follows: • We present a novel two-stage numerical approach for solving stochastic epidemic models, especially SPDEs, that integrates Runge-Kutta schemes with exponential integrators. • We develop a diffusive SIQR model with a general incidence rate and stochastic influences to more precisely represent actual epidemic dynamics. • We prove the mean-square stability and consistency of the suggested numerical scheme to guarantee its dependability in long-term simulations. • We reduce computational complexity and improve efficiency by developing a strategy focusing on time-dependent terms and avoiding spatial discretization. • We solve both deterministic and stochastic versions of the SIQR model by the Euler- Maruyama method, the Nonstandard Finite Difference (NSFD) method, and the proposed scheme, and demonstrate that the proposed method yields lower error and higher accuracy. Infectious diseases have seriously threatened humankind for a long time, and they continue to do so now. The World Health Organization reports that between 1 January 2020 and 1 May 2020, infectious diseases claimed the lives of more than 4 million people across the globe. Consequently, it is crucial to regulate the spread of diseases and have a firm grasp of epidemiological patterns. In recent years, mathematical epidemic modelling has emerged as a powerful and significant tool for comprehending the spread of contagious illnesses. Mathematical models M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 3 of 21 of epidemiology, initially proposed by Kermack and McKendrick, describe the processes that cause epidemics to spread. Both theoretically and pragmatically, their work serves as a reference for future mathematical modelling studies and applications in epidemiology. A plethora of subsequent epidemic models have been put forth in the literature [1–4]. Vaccination, transient immunity, and fluctuating population sizes were all factors in the SIS epidemic model that Li and Ma examined [5]. Kilicman and Hamdan suggested a fractional-order SIR epidemic dengue transmission model [6]. Xiang et al. demonstrated local stability in their discrete SIRS epidemic model with vaccination [7]. Quarantine was thought to be the most essential way to prevent the spread of conta- gious diseases for a long time. Coronavirus, influenza, rubella, and other epidemics have been contained with its help. Quarantine has recently proven effective in limiting the spread of the deadly new coronavirus [8, 9]. As a result, research into how quarantine affects epidemic behaviour is crucial. To illustrate this, Heathcote et al. created an SIQR model [10]. The effects of quarantine on influenza transmission dynamics were investi- gated by Erdem et al. [11]. An SIQR epidemic model was created by Ma et al. [12] that integrated hybrid techniques for vaccination and elimination. An infectious disease’s incidence rate is defined as the number of newly reported cases in a population as a function of time. It is crucial to model the spread of infectious diseases [13–15]. The standard incidence rate βSI and bilinear incidence rate βSI are commonly used in epidemiological models. However, these may be inadequate in scenarios such as population saturation [16, 17]. Non-linear incidence models, such as βIpSq and αI/(S + βI), better capture realistic infection patterns [18–20]. Infectious diseases inflict not only physical suffering on humans but also result in property loss and significant environmental impact. With globalization and technological advancement, cross-regional disease transmission is accelerating. Though many studies on conventional compartment models assume homogeneous mixing [21–25], these assumptions fail in light of heterogeneous disease patterns seen in SARS and AIDS [26, 27]. Many real-world systems, such as the World Wide Web [28, 29] and human contact networks, are scale-free, as described by Barabási and Albert [30]. The scale-free network model leads to power-law degree distributions and challenges traditional threshold-based epidemic predictions [31]. For example, Pastor-Satorras and Vespignani showed that epi- demic thresholds tend to zero in large scale-free networks [32, 33]. Li et al. analyzed a SIQRS model on such networks and highlighted the combined role of quarantine and net- work topology [34]. Demographics and vaccination further influence dynamics, as shown by Huang et al. [35]. A stochastic diffusive model is created in this study [36] to handle the intricate dynamics of CO2 concentration, population expansion, and energy produc- tion. The researcher offers a new computational technique for solving deterministic and stochastic PDEs in the papers [37, 38]. To account for individual behavior changes, Xiao and Ruan proposed a non-monotone incidence rate g(I)S = kIS 1+αI2 , further studied by Li [39]. Random environmental ef- fects, modeled as white noise, influence disease dynamics, prompting the development of stochastic models. These models explore extinction, persistence, and ergodic stationary distributions [40–46]. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 4 of 21 This work offers a new computational method tailored for the stochastic SIQR model with generalized incidence rate. The proposed scheme preserves positivity and bounded- ness—key for epidemiological validity. The main goals of this paper are: • To build and apply a robust numerical method for the stochastic SIQR model. • To study how stochasticity and varying incidence rates affect disease dynamics. • To validate theoretical results with numerical experiments. We solve the SIQR model using four techniques, including the pdepe solver in MAT- LAB, suitable for diffusive and reactive models. As no exact solution exists, the pdepe out- put serves as the reference. The comparison indicates that our proposed scheme achieves higher accuracy than the NSFD method. The remainder of the paper is organized as follows: Section 2 constructs the numer- ical scheme for stochastic differential equations. Section 3 presents a stability analysis. Section 4 discusses the stochastic SIQR model and its local stability. Section 5 provides numerical results. Section 6 concludes the paper. 2. Proposed Numerical Scheme The proposed scheme is explicit and constructed on two time levels. The solution is found by using a predictor-corrector strategy. The predictor stage finds a solution at an arbitrary time level, and the solution obtained by the predictor is utilized in the corrector stage to find the solution at the next time level. The scheme will be proposed for stochastic differential equations but is first developed for deterministic models. To propose the scheme, consider the differential equation: ∂p ∂t = γ ∂2p ∂x2 + f(p) (1) where γ is a constant diffusion coefficient and f(p) is a non-linear source term. 2.1. Predictor-Corrector Framework We propose a two-stage predictor-corrector method in time: Predictor Stage: P̄n+1 i = 1 2 Pn i e ∆t + 1 2 Pn i e 5∆t + (e∆t + e5∆t − 2) 6 ( ∂P ∂t ∣∣∣∣n i − 3Pn i ) (2) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 5 of 21 Corrector Stage: Pn+1 i = aPn i + bP̄n+1 i + c(e∆t − 1) ( ∂P̄ ∂t ∣∣∣∣n+1 i ) (3) Here, a, b, and c are parameters determined using Taylor expansion for improved accuracy. 2.2. Determine Parameters via Taylor Expansion We use the Taylor series expansion: Pn+1 i = Pn i +∆t ∂P ∂t ∣∣∣∣n i + ∆t2 2 ∂2P ∂t2 ∣∣∣∣n i +O(∆t3) (4) Substitute Eqs. (2) and (4) into (3) to match terms. This leads to the following system of equations: 1 = a+ b ( 1 2 e∆t + 1 2 e5∆t − (e∆t + e5∆t − 2) 2 ) ∆t = b · (e ∆t + e5∆t − 2) 6 + c(e∆t − 1) ( 1 2 e∆t + 1 2 e5∆t − (e∆t + e5∆t − 2) 2 ) ∆t2 2 = c(e∆t − 1) · (e ∆t + e5∆t − 2) 6 Solving the system yields the values of a, b, and c. 2.3. Final Deterministic Scheme After computing a, b, and c, the deterministic scheme becomes: Final Predictor: P̄n+1 i = 1 2 Pn i e ∆t + 1 2 Pn i e 5∆t + (e∆t + e5∆t − 2) 6 ( γ ∂2P ∂x2 ∣∣∣∣n i + f(Pn i )− 3Pn i ) (5) Final Corrector: Pn+1 i = aPn i + bP̄n+1 i + c(e∆t − 1) ( γ ∂2P̄ ∂x2 ∣∣∣∣n+1 i + f(P̄n+1 i ) ) (6) 2.4. Space Discretization Using Compact Scheme The above scheme discretizes only in time. For spatial discretization, we apply a compact finite difference scheme: M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 6 of 21 Predictor with Space Discretization: P̄n+1 i = 1 2 Pn i e ∆t + 1 2 Pn i e 5∆t + (e∆t + e5∆t − 2) 6 ( γA−1BPn i + f(Pn i )− 3Pn i ) (7) Corrector with Space Discretization: Pn+1 i = aPn i + bP̄n+1 i + c(e∆t − 1) ( γA−1BP̄n i + f(P̄n+1 i ) ) (8) Here, A and B are matrices defined from the following compact scheme: α1P ′′ i−1 + P ′′ i + α1P ′′ i+1 = b0 Pn i+1 − 2Pn i + Pn i−1 ∆x2 + b1 Pn i+2 − 2Pn i + Pn i−2 4∆x2 (9) with b0 = 4 3 (1− α1), b1 = 1 3 (10α1 − 1) 2.5. Extension to Stochastic Model Now we extend the scheme to the stochastic differential equation: ∂P = ( γ ∂2P ∂x2 + f(P ) ) dt+ σ1PdW (10) where W is a Wiener process and σ1 is the noise intensity. The predictor stage remains the same as in the deterministic case. The corrector stage becomes: Stochastic Corrector: Pn+1 i = aPn i + bP̄n+1 i + c(e∆t − 1) ( γA−1BP̄n+1 i + f(P̄n+1 i ) ) + σ1P n i ∆W (11) where ∆W ∼ N (0,∆t) is sampled from a normal distribution. Summary: Our scheme is explicit and predictor-corrector-based. It handles nonlinearity and diffusion using exponential weights and accurate coefficient matching. It is constructed first for the deterministic model and then extended to the stochastic model by adding the noise term in the corrector stage. Spatial terms are discretized using a compact finite difference scheme, increasing accuracy without requiring a fine mesh and thus improving computational performance. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 7 of 21 3. Stability Analysis Von Neumann analysis is one of the classical techniques for determining the stability conditions of finite difference schemes applied to linear partial differential equations. This criterion, based on Fourier series, transforms linear difference equations into trigonometric equations, from which stability conditions are derived. To apply this method, we use the following transformations: AeiIθ = α1e (i−1)Iθ + eiIθ + α1e (i+1)Iθ (12) BeiIθ = b0 ∆x2 (e(i+1)Iθ − 2eiIθ + e(i−1)Iθ) + b1 4∆x2 (e(i+2)Iθ − 2eiIθ + e(i−2)Iθ) (13) Substituting (12) and (13) into the predictor stage of the scheme (2) with f = 0, we obtain: P̄n+1 i = 1 2 Pn i e ∆t+ 1 2 Pn i e 5∆t+ (e∆t + e5∆t − 2) 6 { γ [4b0(cos θ − 1) + b1(cos 2θ − 1)] 2∆x2(2α1 cos θ + 1) − 3 } Pn i (14) Let β denote the amplification factor from the predictor: P̄n+1 i = βPn i (15) Then, β = 1 2 e∆t + 1 2 e5∆t + (e∆t + e5∆t − 2) 6 { γ [4b0(cos θ − 1) + b1(cos 2θ − 1)] 2∆x2(2α1 cos θ + 1) − 3 } (16) Now, substituting (12) and (13) into the corrector stage of the scheme (with f = 0), we get: Pn+1 i = aPn i + bP̄n+1 i + c(e∆t − 1) { γ [4b0(cos θ − 1) + b1(cos 2θ − 1)] 2∆x2(2α1 cos θ + 1) } P̄n+1 i + σ1P n i ∆W (17) Using (15), we rewrite the above as: Pn+1 i = [ a+ bβ + c(e∆t − 1) { γ [4b0(cos θ − 1) + b1(cos 2θ − 1)] 2∆x2(2α1 cos θ + 1) } β ] Pn i + σ1P n i ∆W (18) Define: β1 = a+ bβ + c(e∆t − 1) { γ [4b0(cos θ − 1) + b1(cos 2θ − 1)] 2∆x2(2α1 cos θ + 1) } β (19) Then, the amplification factor becomes: Pn+1 i Pn i = β1 + σ1∆W (20) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 8 of 21 Taking expectations on the squared modulus yields: E ∣∣∣∣Pn+1 i Pn i ∣∣∣∣2 ≤ 2|β1|2 + 2σ2 1E|∆W |2 (21) If |β1|2 < 1 2 and 2σ2 1 = λ, we obtain: E ∣∣∣∣Pn+1 i Pn i ∣∣∣∣2 ≤ 1 + λ∆t (22) Hence, the proposed scheme is conditionally stable in the mean-square sense. Theorem 1. The proposed schemes (2) and (3) with compact spatial discretization are consistent in the mean-square sense. Proof. Let H be a smooth function. Define: L(H)ni = H((n+ 1)∆t, i∆x)−H(n∆t, i∆x)− γ ∫ (n+1)∆t n∆t Hxx(s, i∆x)ds− σ1 ∫ (n+1)∆t n∆t H(s, i∆x)dW (s) (23) Ln i H = H((n+ 1)∆t, i∆x)−H(n∆t, i∆x)− γb (e∆t + e5∆t − 2) 6 A−1BH(n∆t, i∆x) − γc(e∆t − 1)A−1BH̄((n+ 1)∆t, i∆x)− σ1H(n∆t, i∆x)(W ((n+ 1)∆t)−W (n∆t)) (24) where H̄((n+1)∆t, i∆x) = 1 2 H((n+1)∆t, i∆x) e∆t + 1 2 H(n∆t, i∆x) e5∆t + e∆t + e5∆t − 2 6 { A−1BH(n∆t, i∆x)− 1 2 H(n∆t, i∆x) } (25) Then the error between exact and numerical evolution becomes: E|L(H)ni − Ln i H|2 ≤ 2γ2E ∣∣∣∣∣ ∫ (n+1)∆t n∆t Hxx(s, i∆x)ds− b (e∆t + e5∆t − 2) 6 A−1BH(n∆t, i∆x) +c(e∆t − 1)A−1BH̄((n+ 1)∆t, i∆x) ∣∣2 + 2σ2 1∆t ∫ (n+1)∆t n∆t E [ |H(s, i∆x)−H(n∆t, i∆x)|2 ] ds (26) Using stochastic integral inequalities such as: E ∣∣∣∣∫ t t0 u(s)dW (s) ∣∣∣∣2m ≤ (t− t0) m−1[m(2m− 1)]m ∫ t t0 E|u(s)|2mds (27) Taking the limit as ∆x → 0, ∆t → 0, and (n∆t, i∆x) → (t, x), we get: E|L(H)ni − Ln i H|2 → 0 (28) Thus, the proposed scheme is consistent in the mean-square sense. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 9 of 21 4. Mathematical Model The population is divided into four compartments: susceptible (S(t, x)), infected (I(t, x)), quarantined (Q(t, x)), and recovered (R(t, x)), where t and x represent time and space, respectively. This formulation captures both temporal and spatial dynamics of disease transmission with quarantine and recovery effects. 4.1. Deterministic Diffusive Model The deterministic reaction-diffusion model is given by: ∂S ∂t = d1 ∂2S ∂x2 + Λ− β1SI 1 + αI − dS − ρS (29) ∂I ∂t = d2 ∂2I ∂x2 + β1SI 1 + αI − (γ1 + δ + d+ α2 + ζ)I (30) ∂Q ∂t = d3 ∂2Q ∂x2 + δI − (d+ α3 + ϵ)Q (31) ∂R ∂t = d4 ∂2R ∂x2 + ϵQ+ ρS − dR (32) Here: • Λ is the recruitment rate. • β1 is the transmission rate. • α, α2, α3 denote saturation and comorbidity-induced death. • d is the natural death rate. • δ, ϵ, γ1, ζ represent quarantine, recovery, and disease-induced death rates. • ρ is the vaccination rate. • di (i = 1, 2, 3, 4) are diffusion coefficients. The boundary conditions are: ∂S ∂x ∣∣∣∣ x=0,L = ∂I ∂x ∣∣∣∣ x=0,L = ∂Q ∂x ∣∣∣∣ x=0,L = ∂R ∂x ∣∣∣∣ x=0,L = 0 (33) These Neumann (no-flux) conditions assume no migration at the domain boundaries. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 10 of 21 4.2. Disease-Free Equilibrium (DFE) The disease-free equilibrium (DFE) is obtained by setting the right-hand sides of Eqs. (29)–(32) to zero and assuming I = Q = 0: Λ− dS − ρS = 0 (34) δI = 0 (35) (d+ α3 + ϵ)Q = 0 (36) ϵQ+ ρS − dR = 0 (37) Solving yields the DFE point: E◦ = ( S = Λ d+ ρ , I = 0, Q = 0, R = ρΛ d(d+ ρ) ) (38) 4.3. Local Stability of the DFE Theorem 2. The DFE is locally asymptotically stable if the following condition is satisfied: β1Λ < ζα2d+ d2 + dδ + dγ1 − ζα2ρ+ dρ+ δρ+ γ1ρ Proof. With di = 0, the Jacobian matrix J of the system (29)–(32) is: J =  −d− β1I (1+αI)2 − ρ − β1S (1+αI)2 + 2αβ1IS (1+αI)3 0 0 β1I (1+αI)2 −ζα2 − d− δ − γ1 + β1S (1+αI)2 − 2αβ1IS (1+αI)3 0 0 0 δ −α3 − d− ϵ 0 ρ γ1 ϵ −d  Evaluated at E◦: J |E◦ =  −d− ρ − β1Λ d+ρ 0 0 0 −ζα2 − d− δ + β1Λ d+ρ 0 0 0 δ −α3 − d− ϵ 0 ρ γ1 ϵ −d  The eigenvalues of J |E◦ are: λ1 = −d, λ2 = −α3 − d− ϵ, λ3 = −d− ρ, λ4 = 1 d+ ρ [β1Λ− (stability threshold)] If all eigenvalues are negative, the DFE is locally asymptotically stable. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 11 of 21 4.4. Stochastic Diffusive Model To capture random fluctuations, we introduce stochasticity into the model: dS = ( d1 ∂2S ∂x2 + Λ− β1SI 1 + αI − dS − ρS ) dt+ σ1SdW (39) dI = ( d2 ∂2I ∂x2 + β1SI 1 + αI − (γ1 + δ + d+ α2 + ζ)I ) dt+ σ2IdW (40) dQ = ( d3 ∂2Q ∂x2 + δI − (d+ α3 + ϵ)Q ) dt+ σ3QdW (41) dR = ( d4 ∂2R ∂x2 + ϵQ+ ρS − dR ) dt+ σ4RdW (42) Here, σi (i = 1, 2, 3, 4) are noise intensities, and W represents a Wiener process (Brow- nian motion). 4.5. Model Relevance to Real Epidemics This model integrates several realistic epidemic mechanisms: • Spatial diffusion accounts for movement/contact patterns, critical for modeling localized outbreaks. • Nonlinear incidence rate reflects saturation effects when infection levels are high. • Stochastic terms introduce variability and environmental uncertainty into epi- demic dynamics. • Quarantine and vaccination model disease control strategies. • Stability analysis provides threshold criteria for disease eradication, assisting in public health decision-making. 5. Results and Discussions A computational framework is presented for solving both deterministic and stochastic models. The technique is comprised of two distinct phases, referred to as the predictor and corrector stages. The predictor stage is utilized to ascertain solutions at any given time level. However, the corrector stage determines solutions at the subsequent time level. The predictor step disregards the stochastic component of stochastic partial differential equations, whereas the corrector stage incorporates this stochastic element. The approach provides a second-order accurate solution for deterministic models. The initial stage of the scheme diverges from the conventional Runge-Kutta method, but the subsequent stage aligns with the second stage of the Runge-Kutta scheme; hence, the entire scheme repre- sents a hybrid of a modified exponential integrator and the Runge-Kutta method. The M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 12 of 21 suggested technique is preferable to the existing Euler-Maruyama method due to its supe- rior accuracy for deterministic models while addressing the stochastic component, similar to the current Euler-Maruyama method. For simulation, we take the birth/death rate d = 0.1, transmission rate β1 = 0.1, recruitment rate Λ = 0.1, vaccination rate ρ = 0.1, quarantine rate δ = 0.1, nonlinear incidence factor α = 0.3, comorbidity death from infection α2 = 0.1, comorbidity death from quarantine α3 = 0.1, total population size N = 15 obtained using the initial dis- tribution as S0 + I0 + Q0 + R0 = 3 + 4 + 3 + 5, disease elimination from infected class ζ = 0.1, recovery rate from infection γ1 = 0.3, and removal rate from quarantine ϵ = 0.1. Furthermore, σi = 0.1, di = 0.1 for i = 1, 2, 3, 4. To determine the eigenvalues of the Jacobian matrix at the disease-free equilibrium (DFE) E◦ for the given parameter values, note that at the equilibrium state, S = Λ/(d+ρ) and I = 0, Q = 0, R = ρΛ/(d(d+ ρ)), which gives: J |E◦ =  −0.2 −0.05 0 0 0 −0.16 0 0 0 0.1 −0.3 0 0.1 0.3 0.1 −0.1  Solving the characteristic equation |J − λI| = 0, i.e.,∣∣∣∣∣∣∣∣ −0.2− λ −0.05 0 0 0 −0.16− λ 0 0 0 0.1 −0.3− λ 0 0.1 0.3 0.1 −0.1− λ ∣∣∣∣∣∣∣∣ = 0 we obtain the eigenvalues: λ1 = −0.1, λ2 = −0.2, λ3 = −0.3, λ4 = −0.16 Since all eigenvalues are real and negative, the disease-free equilibrium is locally asymp- totically stable. This supports the theoretical result stated in Theorem 2 for the given parameter set. Comparison of Numerical Schemes (Proposed vs Euler-Maruyama) To compare the performance of the proposed scheme for the stochastic SIQR model, Figure 1 shows the epidemic dynamics obtained using the Euler-Maruyama method (right) and the proposed numerical scheme (left) for the given parameter set: d1 = d2 = d3 = d4 = 0.1, β1 = 0.1, ρ = 0.5, δ = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, x = 0.1429, σ1 = σ2 = σ3 = σ4 = 0.1. The difference in the solutions obtained using the two methods is due to their dif- fering order of accuracy for solving the deterministic part of the model. As previously discussed, the Euler-Maruyama method is first-order accurate, while the proposed scheme is second-order accurate in time. Both methods employ the same approach to incorporate M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 13 of 21 Figure 1: Comparison of proposed and Euler-Maruyama method for solving the stochastic model using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, β1 = 0.1, ρ = 0.5, δ = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, x = 0.1429, σ1 = 0.1, σ2 = 0.1, σ3 = 0.1, σ4 = 0.1. the stochastic (Wiener) terms. The Euler-Maruyama method results in visibly noisier tra- jectories, particularly for the susceptible (purple) and recovered (cyan) subpopulations. In contrast, the proposed scheme yields smoother and more stable dynamics, following ex- pected epidemic trends closely. This highlights the superiority of the proposed method for accurately and efficiently solving complex stochastic partial differential equation models in disease transmission analysis. Effect of Transmission Rate β1 on Susceptible and Infected Classes Next, we investigate the effect of the transmission rate β1 on the susceptible S and infected I classes. To see the effect more clearly, we demonstrate using the deterministic model, i.e., σi = 0, (i = 1, 2, 3, 4). The left subplot shows the susceptible S population over time for different values of β1. The right subplot shows the corresponding infective I population over time. Three values of transmission rate β1 are tested: • Solid lines: β1 = 0.1 • Dashed lines: β1 = 0.3 • Dotted lines: β1 = 0.5 The parameter settings are d1 = d2 = d3 = d4 = 0.1, ρ = 0.3, δ = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1. As β1 increases, the susceptible population declines more rapidly. Higher transmission rates cause faster infection to spread, meaning more people move from susceptible to infective in a shorter time. As β1 increases, the infective population rises sharply at the beginning, peaks earlier, and then declines. Interestingly, even with higher β1, the infective population eventually declines to near zero, indicating that the model includes effective recovery or control mechanisms. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 14 of 21 Figure 2: Effect of transmission rate (from susceptible to infective) on susceptible and infective using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, ρ = 0.3, δ = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. Effect of Quarantine Rate δ on Infective and Quarantined Classes Figure 3 illustrates the effect of the quarantine rate of infective individuals δ on the infective and quarantined compartments in the deterministic SIQR model. The left plot shows the time evolution of the infective I population. The right plot shows the time evolution of the quarantined Q population. Three values of δ are tested: • Solid line: δ = 0.1 • Dashed line: δ = 0.3 • Dotted line: δ = 0.5 The parameter settings are d1 = d2 = d3 = d4 = 0.1, ρ = 0.3, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1. Figure 3: Effect of rate of quarantine of infective on infective and quarantined individuals using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, ρ = 0.3, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. As δ increases (stronger quarantine enforcement), the infective population declines more rapidly. Higher quarantine rates move individuals more quickly out of the infective M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 15 of 21 class into the quarantine class. Increasing δ leads to a higher and earlier peak in the quar- antined population. Over time, the quarantined population also declines due to recovery ϵ or death (α3, d). Influence of Vaccination Rate ρ on Susceptible and Recovered Classes Figure 4 illustrates the influence of the vaccination rate ρ on the susceptible S and recovered R populations in the deterministic SIQR model. The left plot shows how the susceptible population S(t) evolves over time. The right plot shows how the recovered population R(t) evolves over time. Three values of ρ are tested: • Solid line: ρ = 0.1 • Dashed line: ρ = 0.3 • Dotted line: ρ = 0.5 The parameter settings are d1 = d2 = d3 = d4 = 0.1, δ = 0.1, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1. Figure 4: Effect of vaccination rate on susceptible and recovered individuals using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, δ = 0.1, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. As ρ increases, the susceptible population decreases faster. Higher vaccination rates remove individuals from the susceptible class and transfer them to the recovered class. Thus, higher ρ leads to a quicker and higher peak in the recovered population. The recovered curve rises steeply and reaches a higher maximum when vaccination is more aggressive. Comparison with the Nonstandard Finite Difference (NSFD) Method Finally, we compare the performance of the proposed numerical scheme with the non- standard finite difference (NSFD) method. For clarity, in Figures 5 to 8, we only show the results for the deterministic model. The superiority of the proposed scheme for the SDE model is demonstrated in Figure 1 above. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 16 of 21 Figure 5: Comparison of proposed and existing NSFD method for susceptible using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, δ = 0.1, ρ = 0.5, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. Since the exact solution is not available, we use the numerical solution obtained using the MATLAB built-in utility (which can solve parabolic PDEs) as a surrogate for the exact solution. Using d1 = d2 = d3 = d4 = 0.1, δ = 0.1, ρ = 0.5, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, Figures 5–8 show the error incurred by the proposed and NSFD schemes. Figure 6: Comparison of proposed and existing NSFD method for infective using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, δ = 0.1, ρ = 0.5, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. The proposed scheme provides a considerably better approximation than the NSFD method, which suffers from a consistency issue, as discussed in [47]. Moreover, the NSFD method does not allow a higher-order scheme for space discretization. The proposed scheme is not subject to the same limitation. 6. Conclusion In this paper, we suggested a new two-stage computing methodology that combines Runge-Kutta methods with exponential integrator techniques to solve stochastic diffu- sive SIQR epidemic models with a broad incidence rate. The fundamental benefit of the suggested approach is its capacity to discretize just the time-dependent parts of the un- derlying stochastic partial differential equations, lowering computational complexity while M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 17 of 21 Figure 7: Comparison of proposed and existing NSFD method for quarantined individuals using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, δ = 0.1, ρ = 0.5, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. Figure 8: Comparison of proposed and existing NSFD method for recovered using d1 = 0.1, d2 = 0.1, d3 = 0.1, d4 = 0.1, δ = 0.1, ρ = 0.5, β1 = 0.1, γ1 = 0.1, d = 0.1, Λ = 0.1, α = 0.3, α2 = 0.1, α3 = 0.1, ϵ = 0.1, ζ = 0.1, σ1 = 0, σ2 = 0, σ3 = 0, σ4 = 0. preserving great accuracy. A mixture of two schemes has been proposed to handle deterministic and stochastic partial differential equations. The scheme was comprised of two stages called the predictor and corrector stages. The predictor stage found the solution at an arbitrary time level, and the corrector stage found the solution at the next time level. The scheme has been applied to the deterministic and stochastic SIQR model. In addition to this, the deterministic model was solved using the MATLAB solver pdepe. The solver can be used to find solutions for diffusive epidemic models using built-in code available in MATLAB software. We established the stability and consistency of the proposed scheme in the mean-square sense, confirming its theoretical soundness. Numerical experiments were conducted to compare the performance of the proposed method with the widely used Euler-Maruyama method and the nonstandard finite difference (NSFD) scheme. The results demonstrated that the proposed method outperforms existing techniques by producing more accurate solutions with lower numerical errors. The following points can be concluded: • The proposed scheme performed better than the existing nonstandard finite differ- ence scheme in accuracy. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 18 of 21 • The rise in transmission parameters from susceptible to infective produced a decline in susceptible and growth in infective. • The vaccination rate yields a decline in susceptible and growth in recovered individ- uals. This study offers a powerful and effective computational instrument for stochastic epi- demic models characterized by spatial diffusion and non-linear dynamics. The suggested framework can be adapted to many stochastic models in epidemiology and related disci- plines, presenting a viable avenue for future study in computational disease modelling. Acknowledgements The authors would like to acknowledge the support of Prince Sultan University for paying the Article Processing Charges (APC) of this publication. References [1] J. Mann and M. Roberts. Modelling the epidemiology of hepatitis b in new zealand. Journal of Theoretical Biology, 269(1):266–272, 2011. [2] C. Tian, Q. Zhang, and L. Zhang. Global stability in a networked sir epidemic model. Applied Mathematics Letters, 107:106444, 2020. [3] X. Meng, S. Zhao, T. Feng, and T. Zhang. Dynamics of a novel non-linear stochas- tic sis epidemic model with double epidemic hypothesis. Journal of Mathematical Analysis and Applications, 433(1):227–242, 2016. [4] A. El Koufi, J. Adnani, A. Bennar, and N. Yousfi. Analysis of a stochastic sir model with vaccination and non-linear incidence rate. International Journal of Differential Equations, 2019, 2019. [5] J. Li and Z. Ma. Qualitative analyses of sis epidemic model with vaccination and vary- ing total population size. Mathematical and Computer Modelling, 35(11–12):1235– 1243, 2002. [6] A. Kilicman. A fractional order sir epidemic model for dengue transmission. Chaos, Solitons & Fractals, 114:55–62, 2018. [7] L. Xiang, Y. Zhang, and J. Huang. Stability analysis of a discrete sirs epidemic model with vaccination. Journal of Difference Equations and Applications, 26(3):309–327, 2020. [8] T. Odagaki. Exact properties of siqr model for covid-19. Physica A, 564:125564, 2021. [9] T. Odagaki. Analysis of the outbreak of covid-19 in japan by siqr model. Infectious Disease Modelling, 5:691–698, 2020. [10] H. Hethcote, M. Zhien, and L. Shengbing. Effects of quarantine in six endemic models for infectious diseases. Mathematical Biosciences, 180(1–2):141–160, 2002. [11] M. Erdem, M. Safan, and C. Castillo-Chavez. Mathematical analysis of an siqr influenza model with imperfect quarantine. Bulletin of Mathematical Biology, 79(7):1612–1636, 2017. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 19 of 21 [12] Y. Ma, J.B. Liu, and H. Li. Global dynamics of an siqr model with vaccination and elimination hybrid strategies. Mathematics, 6(12):328, 2018. [13] H. Cheng, F. Wang, and T. Zhang. Multi-state dependent impulsive control for holling i predator-prey model. Discrete Dynamics in Nature and Society, 2012, 2012. [14] W. Zou, J. Xie, and Z. Xiong. Stability and hopf bifurcation for an eco-epidemiology model with holling-iii functional response and delays. International Journal of Biomathematics, 1(3):377–389, 2008. [15] Y. Jin, W. Wang, and S. Xiao. An sirs model with a non-linear incidence rate. Chaos, Solitons & Fractals, 34(5):1482–1497, 2007. [16] M.A. Khan, Y. Khan, and S. Islam. Complex dynamics of an seir epidemic model with saturated incidence rate and treatment. Physica A, 493:210–227, 2018. [17] E.A. Algehyne and R. ud Din. On global dynamics of covid-19 by using sqir type model under non-linear saturated incidence rate. Alexandria Engineering Journal, 60(1):393–399, 2021. [18] S. Han and C. Lei. Global stability of equilibria of a diffusive seir epidemic model with non-linear incidence. Applied Mathematics Letters, 98:114–120, 2019. [19] A. Kumar and M. Kumar. A study on the stability behavior of an epidemic model with ratio-dependent incidence and saturated treatment. Theory in Biosciences, 139(2):225–234, 2020. [20] Y. Zhang, K. Fan, S. Gao, Y. Liu, and S. Chen. Ergodic stationary distribution of a stochastic sirs epidemic model incorporating media coverage and saturated incidence rate. Physica A, 514:671–685, 2019. [21] H. Hethcote. The mathematics of infectious diseases. SIAM Review, 42:599–653, 2000. [22] G. Zhu, G. Chen, H. Zhang, and X. Fu. Propagation dynamics of an epidemic model with infective media connecting two separated networks of populations. Communi- cations in Nonlinear Science and Numerical Simulation, 20:240–249, 2015. [23] S. Djilali, B. Ghanbari, S. Bentout, and A. Mezouaghi. Turing-hopf bifurcation in a diffusive mussel-algae model with time-fractional-order derivative. Chaos, Solitons & Fractals, 138:109954, 2020. [24] B. Pgda, D. Cckc, and F. Rmse. Determining the optimal strategy for reopening schools, the impact of test and trace interventions, and the risk of occurrence of a second covid-19 epidemic wave in the uk: a modelling study. The Lancet Child & Adolescent Health, 4:817–827, 2020. [25] A. Alkhazzan, J. Wang, Y. Nie, H. Khan, and J. Alzabut. A stochastic susceptible vaccinees infected recovered epidemic model with three types of noises. Mathematical Methods in the Applied Sciences, 47(11):8748–8770, 2024. [26] F. Liljeros, C.R. Edling, and L.A.N. Amaral. The web of human sexual contacts. Nature, 411:907–908, 2001. [27] Z. Wang, M. Andrews, Z. Wu, L. Wang, and C. Bauch. Coupled disease–behavior dynamics on complex networks: a review. Physics of Life Reviews, 15:1–29, 2015. [28] R. Albert, H. Jeong, and A.L. Barabási. Diameter of the world wide web. Nature, 401:130–131, 1999. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 20 of 21 [29] M. Faloutsos, P. Faloutsos, and C. Faloutsos. On power-law relationships of the internet topology. Computer Communication Review, 29:251–262, 1999. [30] A.L. Barabási and R. Albert. Emergence of scaling in random networks. Science, 286:509–512, 1999. [31] B. Geraldine and M. Antoine. Social interactions and the prophylaxis of si epidemics on networks. Journal of Mathematical Economics, 93:102486, 2021. [32] P.S. Romualdo and V. Alessandro. Epidemic spreading in scale-free networks. Physical Review Letters, 86:3200–3203, 2000. [33] X. Wei, G. Xu, and L. Liu. Global stability of endemic equilibrium of an epidemic model with birth and death on complex networks. Physica A, 477:8–84, 2017. [34] T. Li, Y. Wang, and Z. Guan. Spreading dynamics of a siqrs epidemic model on scale-free networks. Communications in Nonlinear Science and Numerical Simulation, 19:686–692, 2014. [35] S. Huang, F. Chen, and L. Chen. Global dynamics of a network-based siqrs epidemic model with demographics and vaccination. Communications in Nonlinear Science and Numerical Simulation, 43:296–310, 2017. [36] M. S. Arif, K. Abodayeh, H. M. Al-Khawar, and Y. Nawaz. Stochastic diffusive modeling of co2 emissions with population and energy dynamics. Emerging Science Journal, 9(1):210–228, 2025. [37] M. S. Arif. A novel explicit scheme for stochastic diffusive sis models with treatment effects. Partial Differential Equations in Applied Mathematics, page 101215, 2025. [38] M. S. Arif. Numerical analysis of deterministic and stochastic model of covid-19 co-infection with influenza. European Journal of Pure and Applied Mathematics, 18(2):6005–6005, 2025. [39] H. Tahir, A. Din, K. Shah, B. Abdalla, and T. Abdeljawad. Advances in stochastic epidemic modeling: tackling worm transmission in wireless sensor networks. Mathe- matical and Computer Modelling of Dynamical Systems, 30(1):658–682, 2024. [40] M. Mediani, A. Slama, A. Boudaoui, and T. Abdeljawad. Analysis of a stochastic seiuirr epidemic model incorporating the ornstein-uhlenbeck process. Heliyon, 10(16), 2024. [41] Q. Liu, D. Jiang, and N. Shi. Stationary distribution and extinction of a stochastic seir epidemic model with standard incidence. Physica A, 469:510–517, 2017. [42] Y. Zhang and Y. Li. Evolutionary dynamics of stochastic seir models with migration and human awareness in complex networks. Complexity, 2020:3768083, 2020. [43] R. Zhao, Q. Liu, and M. Sun. Dynamical behavior of a stochastic siqs epidemic model on scale-free networks. Journal of Applied Mathematics and Computing, 9:1–26, 2021. [44] Y. Chen, B. Wen, and Z. Teng. The global dynamics for a stochastic sis epidemic model with isolation. Physica A, 492:1604–1624, 2018. [45] L.J. Pan, J.D. Cao, and A. Ahmed. Stability of reaction–diffusion systems with stochastic switching. Nonlinear Analysis: Modelling and Control, 24:315–331, 2019. [46] M.E. Fatini, R. Pettersson, and I. Sekkak. A stochastic analysis for a triple delayed siqr epidemic model with vaccination and elimination strategies. Journal of Applied Mathematics and Computing, 64:781–805, 2020. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6176 21 of 21 [47] Syed Ahmed Pasha, Yasir Nawaz, and Muhammad Shoaib Arif. On the nonstandard finite difference method for reaction–diffusion models. Chaos, Solitons & Fractals, 166:112929, 2023.