EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6292 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Hybrid Finite Difference Approach for Solving Fuzzy Stochastic SIR-β Model with Diffusion and Incidence Rate 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. A finite difference scheme is proposed for solving stochastic fuzzy partial differential equations. The scheme is explicit and constructed on two-time levels. The first stage of the scheme is the modified time integrator. The stability and consistency of the scheme in the mean square sense are also provided. The scheme is applied to the mathematical model of stochastic fuzzy SIR-β model using incidence rate. The scheme is compared with existing nonstandard finite difference schemes for solving deterministic models. The scheme performs better than the existing nonstandard finite difference method in most compared figures. Using a computationally efficient and robust system, this work develops numerical methods for solving epidemiological models under uncertainty, enabling researchers and lawmakers to improve knowledge and control of infectious diseases. 2020 Mathematics Subject Classifications: 35R60, 65M06, 92D30 Key Words and Phrases: Finite Difference Scheme, Exponential Integrator, Stability in Mean Square Sense, SIR-β Model, Incidence Rate, Diffusion in Disease Modeling 1. Introduction Understanding the dynamics of infectious diseases depends much on mathematical modelling, which also helps with control plans and decision-making. Described as the transmission of infectious diseases within a community, the Susceptible-Infected-Recovered (SIR) model is one of the extensively investigated models in epidemiology. Several expan- sions of the SIR model are proposed to grasp more complex real-world scenarios involving stochastic effects, fuzzy uncertainty, diffusion processes, and non-linear incidence rates. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6292 Email addresses: marif@psu.edu.sa (M. S. Arif), kamal@psu.edu.sa (K. Abodayeah), 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), 6292 2 of 24 These developments allow the modelling of uncertainty arising from different population distributions, erroneous data, and environmental fluctuations. In this work, we propose a computational method for a fuzzy stochastic SIR-β model with incidence rate to consider the intrinsic uncertainty in disease transmission. While the stochastic character of the model fits erratic changes in infection and recovery rates, including fuzzy logic lets the model manage imprecise or vague elements like population interactions and disease transmission rates. Including non-linear effects lets the parameter reflect a modified incidence rate, allowing a more flexible approach to disease transmission modelling. This research aims to: • Develop a fuzzy stochastic SIR-β model incorporating random perturbations and imprecise parameters. • Design a computational scheme to solve the proposed model efficiently. • Analyze the effects of stochasticity, fuzziness, and varying incidence rates on disease spread. • Validate the model using numerical simulations and real-world epidemic data. Mathematical models became essential after McKendrick and Kermack [1] ground- breaking work. This study mostly aims to investigate the causes and control of infectious diseases. The erratic changes in the external environment in ecological systems always affect the power of humans. Models of the rate of disease transmission during epidemics have been modified to consider numerous environmental aspects, including the amount of UV light and the fact that many viruses and bacteria live and spread under moist conditions [2–4]. Changes in the range of food sources are influenced by unpredictable oscillations in precipitation, which in turn are influenced by environmental impacts and population growth [5]. One widely used approach in the literature involves randomly switching modes in a finite state such that the continuous-time values of the Markov chain are separated from the values in space. This separation modifies the fundamental parame- ters connected with the Markov chain state flipping [6–8]. A mathematical equation helps the dynamical system to become a piecewise deterministic Markov process. Models of epidemics with deterministic main parameters are unlikely to be applicable to the actual world because they are based on biological systems and subject to stochastic environmen- tal conditions. Analysis of the disease-causing effects of abnormal mutations in naturally occurring mental states on infection transmission in the target population is crucial. De- spite its obvious importance, environmental noise has received scant attention in infectious illness literature. An epidemic SIS model for Markovian switching was first developed by Gray, Greenhalgh, Mao et al. [9]. Some writers have shown that using irregular switch- ing of environmental regimes yields trustworthy findings for time discretization [10, 11]. Many writers have utilized LSCM to approximate solutions of stochastic differential mod- els [12, 13]. The numerical approximation of integrals and differential equations are two other applications of Legendre polynomials [14–18]. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 3 of 24 According to [19], a deterministic SIRS model examines the impact of Markov switch- ing. Because of the COVID-19 pandemic, SIRS stochastic models are now widely used to research the spread of the virus. An improved model for COVID-19 transmission was de- veloped, for instance, by applying the traditional deterministic SIR epidemic model [20]. Using Saudi Arabia as a case study, Khan et al. [21] utilized differential equations of fractional order to model COVID-19 mathematically and offered a complete analysis. Atangana, Seda, and Prolonger’s method was based on integral and piecewise differen- tial operators [22, 23]. The dynamics of COVID-19 are analysed via the fuzzy Caputo approach in [24]. The authors of these works developed mathematical models for COVID- 19 transmission via deterministic methods. System parameters are variable, oscillating around a mean value, which constrains the capacity of these models to represent reality compared to stochastic models accurately [25]. The integrability of a SIRS stochastic differential equation is demonstrated to be as- sociated with its Lie symmetry in [26]. The Lie point symmetries of the stochastic SIRS model are identified using the derived equation. A unique viable solution is demonstrated and the corresponding message-passing system of equations is presented. In the case of a symmetric graph (network) with a homogeneous special condition, a generalized upper limit is established for the projected epidemic size at a specific time t > 0 [27]. In the actual world, exogenous volatility is one of the main factors in the spread of epidemic dis- eases. Since it is crucial to incorporate unpredictability into the model, stochastic SIRS models are seen as more realistic and practical for studying and simulating the evaluation of epidemic diseases. The use of stochastic SIRS models in disease surveillance has been the subject of extensive research [28–32]. The parameters used in epidemic models are not precise enough to account for population- wide differences in vulnerability, exposure, infectiousness, and recovery. Considerations like resistance patterns, age groupings, and demographic types may lead to differing out- comes. More accurate models are needed to account for these different levels of people. According to [33], a new strategy is necessary for epidemic systems, especially those in- volving infectious diseases, because of the high uncertainty. In many fields where human judgment is crucial for evaluation and reasoning, fuzzy sets and fuzzy logic have been widely used to solve real-world problems [34–38]. These domains include, but are not limited to, engineering, economics, medicine, and countless others. Epidemiologists have also made use of this theoretical framework. Barros et al. proposed an SI model [39] that used fuzzy logic and treated the transmission coefficient like a fuzzy set. Ortega et al. [40] employed fuzzy logic for problem prediction concerning infectious illness epidemiology. There was discussion of a rabies model including partially vaccinated dogs. Using fuzzy set theory, Mondal et al. [41] developed a SIS model to analyze the epi- demic. A SIR model was developed and numerically and mathematically tested by Das and Pal [42]. Under fuzzy conditions, Sadhukhan et al. [43] investigated harvesting opti- mization within a food chain model. To better understand the dynamics of coronavirus disease, Li et al. [44] built a fuzzy SEIR model that was strengthened by the confidence in- dex theory. A SIR model with fuzzy parameters was presented by Abdy et al. [45] to show how COVID-19 behaves. Allehiany et al. also explored using Euler, RK-4, and NSFD M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 4 of 24 techniques in a fuzzy SIR model [46]. Mickens [47] first proposed the NSFD methodology, and since then, several scholars have used it to solve systems of differential equations [48– 50]. Using fuzzy logic, Adak and Jana investigated a SIS model that included treatment control. Traditional deterministic epidemic models overlook the inherent uncertainty and ran- domness of disease transmission. Although stochastic SIR models consider random fluc- tuations, estimating their well-defined parameters can be challenging in practical settings. Conversely, epidemic modeling greatly benefits from fuzzy logic’s formal foundation for dealing with unknown or imprecise parameters. An adaptive model that better captures real-world epidemic scenarios can be developed by combining stochasticity and fuzziness. The significance of this research lies in: (i) The fuzzy stochastic method enhances the depiction of uncertainty and better cap- tures the dynamics of actual epidemics by integrating random fluctuations and im- precise data. (ii) In the generalized incidence rate, including the parameter β enables more adaptabil- ity in depicting several disease transmission patterns, including saturation effects and behavioural changes. (iii) Efficient Numerical Solutions for the Complex Fuzzy Stochastic SIR-β Model: This computational scheme allows for practical applications like epidemic forecasting and intervention planning by ensuring efficient numerical solutions. Although traditional deterministic models assume correct input parameters, in real- world epidemiological research, factors like environmental variations, measurement errors, and unpredictability in disease transmission cause uncertainty. While stochastic differ- ential equations (SDEs) provide a structure for including unpredictability in the system, fuzzy logic allows one to characterize imperfect or hazy knowledge. In epidemiological models, stochasticity mixed with fuzziness generates fuzzy stochastic differential equa- tions (FSDEs), more faithfully representing disease dynamics. Numerical methods for handling FSDEs are essential in the lack of analytical solutions for complex models. Standard finite difference methods, nonstandard finite difference schemes, and stochastic numerical techniques are among several computing approaches suggested here. Still, current techniques can suffer from instability, poor precision, and excessive computing expense. Particularly in the framework of an extended SIR-β model with incidence rate and diffusion effects, we offer in this work an effective finite difference scheme for solving stochastic fuzzy partial differential equations. The subsequent sections of this work are structured as follows: Section 2 delineates the mathematical formulation of the fuzzy stochastic SIR-β model. Section 3 describes the construction of the finite difference scheme and its stability analysis. Section 4 provides numerical simulations and comparative results with the existing NSFD method. Finally, Section 5 concludes the study with discussions on future research directions. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 5 of 24 2. Proposed Fuzzy Numerical Scheme 2.1. Fuzzy Set A fuzzy set A in X is a function from set X to the unit interval [0, 1]. A membership function µA(x) characterizes a fuzzy set. Each member in A is an ordered pair consisting of the element in X and its grade of membership. If the membership function is expressed as a triangle, then it is called a triangular membership function. 2.2. Proposed Fuzzy Numerical Scheme We propose a fuzzy numerical scheme to solve partial differential equations with fuzzy parameters. The β1-cut approach is used to represent fuzzy parameters. The time dis- cretization is based on a two-stage integrator, and a compact finite difference method is used for spatial discretization. A fuzzy PDE with a fuzzy parameter γ is written as: ∂v ∂t = γ ∂2v ∂x2 + g(v) (1) where γ is a fuzzy triangular number defined as [0.3β1, 1− 0.5β1]. Applying the β1-cut method: ∂ ∂t [v1(t, x, β1), v2(t, x, β1)] = γ ∂2 ∂x2 [v1(t, x, β1), v2(t, x, β1)] + g ([v1(t, x, β1), v2(t, x, β1)]) (2) 2.3. First Stage of the Scheme Let v = v1. The first stage is defined as: v̄n+1 i = 1 2 vni e 6∆t + 1 2 vni e −5∆t + ( e6∆t + e−5∆t − 2 )( ∂v ∂t ∣∣∣∣n i − 1 2 vni ) (3) 2.4. Second Stage of the Scheme The second stage for the final update is: vn+1 i = avni + bv̄n+1 i + c(e∆t − 1) ∂v̄ ∂t ∣∣∣∣n+1 i (4) Using Taylor expansion: vn+1 i = vni +∆t ∂v ∂t ∣∣∣∣n i + (∆t)2 2 ∂2v ∂t2 ∣∣∣∣n i +O(∆t3) (5) By matching coefficients, we get: 1 = a+ b ( 1 2e 6∆t + 1 2e −5∆t − 1 2(e 6∆t + e−5∆t − 2) ) ∆t = b(e6∆t + e−5∆t − 2) + c(e∆t − 1) ( 1 2e 6∆t + 1 2e −5∆t − 1 2(e 6∆t + e−5∆t − 2) ) (∆t)2 2 = c(e∆t − 1)(e6∆t + e−5∆t − 2) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 6 of 24 2.5. Discretization of Fuzzy Spatial Terms Updating Eqs. (4) and (5) for spatial discretization: v̄n+1 i = 1 2 vni e 6∆t + 1 2 vni e −5∆t + ( e6∆t + e−5∆t − 2 )( 0.3β1 ∂2v ∂x2 ∣∣∣∣n i + g(vni )− 1 2 vni ) (6) vn+1 i = avni + bv̄n+1 i + c(e∆t − 1) ( 0.3β1 ∂2v̄ ∂x2 ∣∣∣∣n+1 i + g(v̄n+1 i ) ) (7) For compact spatial discretization, define matrices A and B: α1v ′′ i−1 + v′′i + α1v ′′ i+1 = a0 vni+1 − 2vni + vni−1 ∆x2 + a1 vni+2 − 2vni + vni−2 4∆x2 (8) with a0 = 4 3 (1− α1), a1 = 1 3 (10α1 − 1) Compact form: v̄n+1 i = 1 2 vni e 6∆t + 1 2 vni e −5∆t + ( e6∆t + e−5∆t − 2 )( 0.3β1A −1Bvni + g(vni )− 1 2 vni ) (9) vn+1 i = avni + bv̄n+1 i + c(e∆t − 1) ( 0.3β1A −1Bv̄n+1 i + g(v̄n+1 i ) ) (10) 2.6. Extension to Fuzzy Stochastic Model We now extend to the stochastic fuzzy model: ∂v = 0.3β1 ∂2v ∂x2 dt+ g(v)dt+ σvdW (11) Here, dW ∼ N (0,∆t) is the Wiener increment. The first stage remains unchanged. The second stage becomes: vn+1 i = avni + bv̄n+1 i + c(e∆t − 1) ( 0.3β1A −1Bv̄n+1 i + g(v̄n+1 i ) ) + σvni ∆W (12) This formulation enhances stability, accuracy, and efficiency in solving fuzzy stochastic partial differential equations. 3. Stability Analysis The Von Neumann stability analysis is a widely used criterion to determine the con- ditions for numerical stability. It applies to both linear and nonlinear partial differential equations (PDEs). For linear PDEs, it provides an exact stability condition, while for non- linear PDEs, it gives an estimate. This analysis involves transforming the finite difference scheme into a trigonometric form and deriving stability constraints. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 7 of 24 The Fourier mode transformations are given as: AeiIθ = α1e i(I−1)θ + eiIθ + α1e i(I+1)θ (13) BeiIθ = a0 ( ei(I−1)θ − 2eiIθ + ei(I+1)θ ) ∆x2 + a1 ( ei(I−2)θ − 2eiIθ + ei(I+2)θ ) 4∆x2 (14) Applying these to the first stage of the scheme with g = 0, we obtain: v̄n+1 i = 1 2 vni e 6∆t + 1 2 vni e −5∆t + (e6∆t + e−5∆t − 2) { 0.3β1 [ 4a0(cos θ − 1) + a1(cos 2θ − 1) 2∆x2(2α1 cos θ + 1) ] − 1 2 } vni (15) Define: δ = 1 2 e6∆t + 1 2 e−5∆t + (e6∆t + e−5∆t − 2) { 0.3β1 [ 4a0(cos θ − 1) + a1(cos 2θ − 1) 2∆x2(2α1 cos θ + 1) ] − 1 2 } (16) So: v̄n+1 i = δvni (17) Substitute into the second stage: vn+1 i = avni + bδvni + c(e∆t − 1) [ 0.3β1 ( 4a0(cos θ − 1) + a1(cos 2θ − 1) 2∆x2(2α1 cos θ + 1) )] δvni + σvni ∆W (18) Let: δ1 = a+ bδ + c(e∆t − 1) [ 0.3β1 ( 4a0(cos θ − 1) + a1(cos 2θ − 1) 2∆x2(2α1 cos θ + 1) )] δ (19) Then the update becomes: vn+1 i = δ1v n i + σvni ∆W (20) and the amplification factor is: vn+1 i vni = δ1 + σ∆W (21) Taking expectations: E ∣∣∣∣vn+1 i vni ∣∣∣∣2 ≤ 2E|δ1|2 + 2σ2E|∆W |2 (22) If |δ1|2 < 1 2 , then: E ∣∣∣∣vn+1 i vni ∣∣∣∣2 ≤ 1 + 2σ2E|∆W |2 (23) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 8 of 24 Let λ = 2σ2, then since E|∆W |2 = ∆t, it follows: E ∣∣∣∣vn+1 i vni ∣∣∣∣2 ≤ 1 + λ∆t (24) Hence, the proposed scheme is conditionally stable in the mean square sense. Theorem 1. The proposed scheme (10) and (14) with g = 0 is consistent in the mean square sense. Proof. Let G be a smooth function. Define: L(G)ni = G((n+ 1)∆t, i∆x)−G(n∆t, i∆x)− 0.3β1 ∫ (n+1)∆t n∆t Gxx(s, i∆x) ds − σ ∫ (n+1)∆t n∆t G(s, i∆x) dW (s) (25) Ln i G = G((n+ 1)∆t, i∆x)−G(n∆t, i∆x)− 0.3β1b(e 6∆t + e−5∆t − 2)A−1BG(n∆t, i∆x) − c 0.3β1(e ∆t − 1)A−1BḠ((n+ 1)∆t, i∆x)− σG(n∆t, i∆x)(W ((n+ 1)∆t)−W (n∆t)) (26) where Ḡ((n+ 1)∆t, i∆x) = 1 2 G((n+ 1)∆t, i∆x)e6∆t + 1 2 G(n∆t, i∆x)e−5∆t + (e6∆t + e−5∆t − 2) [ A−1BG(n∆t, i∆x)− 1 2 G(n∆t, i∆x) ] Subtracting (26) from (25), we have: E |L(G)ni − Ln i G|2 ≤ E ∣∣∣∣∣−0.3β1 ∫ (n+1)∆t n∆t Gxx(s, i∆x) ds− σ ∫ (n+1)∆t n∆t G(s, i∆x) dW (s) +0.3β1b(e 6∆t + e−5∆t − 2)A−1BG(n∆t, i∆x) + 0.3β1(e ∆t − 1)A−1BḠ((n+ 1)∆t, i∆x) +σG(n∆t, i∆x)(W ((n+ 1)∆t)−W (n∆t))|2 (27) Inequality (27) becomes: E |L(G)ni − Ln i G|2 ≤ 2(0.3β1) 2E ∣∣∣∣∣ ∫ (n+1)∆t n∆t Gxx(s, i∆x) ds− b(e6∆t + e−5∆t − 2)A−1BG(n∆t, i∆x) +c(e∆t − 1)A−1BḠ((n+ 1)∆t, i∆x) ∣∣2 + 2σ2E ∣∣∣∣∣ ∫ (n+1)∆t n∆t G(s, i∆x) dW (s)−G(n∆t, i∆x)(W ((n+ 1)∆t)−W (n∆t)) ∣∣∣∣∣ 2 (28) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 9 of 24 Using the Itô isometry inequality: E ∣∣∣∣∫ t t0 u(s, ω) dW (s) ∣∣∣∣2m ≤ (t− t0) m−1[m(2m− 1)]m ∫ t t0 E |u(s, ω)|2m ds (29) Apply inequality (29) to (28) to get: E |L(G)ni − Ln i G|2 ≤ 2(0.3β1) 2E ∣∣∣∣∣ ∫ (n+1)∆t n∆t Gxx(s, i∆x) ds− b(e6∆t + e−5∆t − 2)A−1BG(n∆t, i∆x) +c(e∆t − 1)A−1BḠ((n+ 1)∆t, i∆x) ∣∣2 + 2σ2∆t ∫ (n+1)∆t n∆t E |G(s, i∆x)−G(n∆t, i∆x)|2 ds (30) Taking the limit as ∆x → 0, ∆t → 0, and (n∆t, i∆x) → (t, x), we have: E |L(G)ni − Ln i G|2 → 0 (31) Thus, the proposed scheme is consistent in the mean square sense. 4. Mathematical Model A mathematical model exists for modelling epidemic diseases. Among existing models, the Susceptible-Infective-Recovered (SIR) model is powerful, consisting of susceptible, infective, and recovered categories. This study focuses on the SIR-β model, where the interaction rate β is not constant but evolves dynamically. Effects of diffusion, stochastic perturbations, and fuzzy parameters are incorporated to represent real-world uncertainty in disease transmission. The governing equations incorporating diffusion and incidence rate effects are given as: ∂S ∂t = d1 ∂2S ∂x2 + Λ− βSI 1 +mI − (µ+ δ)S, (32) ∂I ∂t = d2 ∂2I ∂x2 + βSI 1 +mI − (µ+ θ + γ)I, (33) ∂R ∂t = d3 ∂2R ∂x2 + γI + δS − µR, (34) ∂β ∂t = d4 ∂2β ∂x2 + β2 1 + c1I − αβ, (35) subject to no-flux boundary conditions: ∂S ∂x = ∂I ∂x = ∂R ∂x = ∂β ∂x = 0 for x = 0, x = L. (36) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 10 of 24 Here, L is a constant spatial domain length, Λ is the birth rate of new susceptible individuals, µ is the natural death rate, θ is the disease mortality rate, δ represents the immunization effect due to vaccination, and γ is the recovery rate. The incidence rate is given by βSI 1+mI , which includes saturation effects, and the adaptive evolution of the transmission rate β follows the equation β2 1+c1I − αβ. The equilibrium point under the steady-state condition and zero diffusion (d1 = d2 = d3 = d4 = 0) satisfies: Λ− βSI 1 +mI − (µ+ δ)S = 0, (37) βSI 1 +mI − (µ+ θ + γ)I = 0, (38) γI + δS − µR = 0, (39) β2 1 + c1I − αβ = 0. (40) The disease-free equilibrium is: B∗ = ( Λ δ + µ , 0, δΛ µ(δ + µ) , β2 α ) . (41) Theorem 2. The model (32)–(35) with d1 = d2 = d3 = d4 = 0 is locally stable at B∗ if: β2Λ < αδγ + αδµ+ αγµ+ αµ2 + αδθ + αµθ. Proof. The Jacobian matrix of the system (32)–(35) at d1 = d2 = d3 = d4 = 0 is: J =  −δ − βSI 1+mI − µ βmSI (1+mI)2 − βS 1+mI 0 − IS 1+mI βI 1+mI −γ − µ− mβSI (1+mI)2 + βS 1+mI − θ 0 IS 1+mI δ γ −µ 0 0 − β2c1 (1+c1I)2 0 −α  . (42) At equilibrium B∗, the Jacobian reduces to: J |B∗ =  −δ − µ β2Λ α(δ+µ) 0 0 0 −γ − µ+ β2Λ α(δ+µ) − θ 0 0 δ γ −µ 0 0 −β2c1 0 −α  . (43) The eigenvalues are: λ1 = −α, λ2 = −µ, λ3 = −δ−µ, λ4 = −αδγ − αδµ− αγµ− αµ2 + β2Λ− αδθ − αµθ α(δ + µ) . M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 11 of 24 Stability requires λ4 < 0, which gives the condition: β2Λ < αδγ + αδµ+ αγµ+ αµ2 + αδθ + αµθ. This ensures the local asymptotic stability of the disease-free equilibrium. 4.1. Extension to Fuzzy Stochastic Diffusive Model The fuzzy stochastic diffusive model incorporating stochastic effects and fuzzy param- eters is given as follows: dS = ( d1 ∂2S ∂x2 + Λ− βSI 1 +mI − (µ+ δ)S ) dt+ σ1S dW, (44) dI = ( d2 ∂2I ∂x2 + βSI 1 +mI − (µ+ θ + γ)I ) dt+ σ2I dW, (45) dR = ( d3 ∂2R ∂x2 + γI + δS − µR ) dt+ σ3RdW, (46) dβ = ( d4 ∂2β ∂x2 + β2 1 + c1I − αβ ) dt+ σ4β dW, (47) where σ1, σ2, σ3, and σ4 denote intensities of stochastic noise, dW is a Wiener incre- ment, and parameters µ and δ are triangular fuzzy numbers. 5. Results and Discussions We present a detailed simulation study with the following objectives: • The proposed scheme applies to deterministic and stochastic models and is con- structed on two-time levels n and (n+1). The solution at (n+1)th is obtained using the solution at the nth time level and an arbitrary intermediate time level. • First Stage: Computes the solution at the arbitrary time level using the nth time level solution. This stage is identical for both deterministic and stochastic models, does not involve any stochastic variation, and is a modification of the exponential integrator, differing from traditional Runge-Kutta schemes. • Second Stage: Handles the Wiener process term in the stochastic partial differen- tial equation, following the Euler-Maruyama method, which is first-order accurate for deterministic models. However, the proposed scheme improves accuracy by being second-order in time. • When applied with a compact spatial scheme, the scheme is conditionally stable, ensuring better numerical stability and higher accuracy than existing methods. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 12 of 24 5.1. Comparison of Proposed and Euler-Maruyama Method for Stochas- tic Model Figure 1 compares the performance of the proposed numerical scheme and the Euler- Maruyama method for solving the stochastic fuzzy SIR-β model. The comparison is made for a specific set of parameters, which include diffusion coefficients, birth rate, transmission rate, recovery rate, and other epidemiological factors, and their values are defined as follows: d1 = d2 = d3 = d4 = 0.1, µ = 0.1, θ = 0.1, γ = 0.4, c1 = 0.1, Λ = 0.1, β1 = 0.1, α = 0.1, m = 0.1. The left panel represents the results obtained using the proposed numerical scheme. The right panel represents the results obtained using the Euler-Maruyama method. Both plots display the time evolution of the four key variables: • Susceptible S – Blue Curve: Initially, the number of susceptible individuals is high. Over time, as individuals become infected, the susceptible population de- creases. In both methods, the susceptible population drops to near zero. • Infective I – Yellow Curve: Initially low but rapidly increases as the infection spreads. Eventually, it decreases as more individuals recover or die. The two meth- ods exhibit similar infection peaks, but the proposed method shows slightly more stability. • Recovered R – Magenta Curve: The number of recovered individuals increases as infected individuals recover. The proposed method produces a smooth recovery curve compared to the Euler-Maruyama method, which shows more fluctuations due to stochastic noise. • Transmission rate β(t) – Cyan Curve: The transmission rate is dynamic and influenced by stochastic variations. In both methods, β(t) starts high and fluctuates over time due to random variations. The proposed method controls fluctuations better, whereas Euler-Maruyama introduces more noise. The proposed scheme demonstrates a smoother transition and better numerical sta- bility. The Euler-Maruyama method introduces higher stochastic variations, leading to greater fluctuations in the infected and recovered populations. The proposed method is second-order accurate in time, whereas the Euler-Maruyama method is only first-order ac- curate, making the proposed method more precise. The stability of the proposed method ensures that the disease-free equilibrium is better preserved, as seen in the declining sus- ceptible and infected populations. A quantitative comparison is shown in Table 1, comparing the proposed scheme and the Euler-Maruyama method based on key performance metrics. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 13 of 24 Table 1: Quantitative comparison between the proposed scheme and the Euler-Maruyama method. Method Mean Final Peak Time to Mean Final Fluct. in Comp. Susceptible Infection Level Peak Recovered β(t) Time (s) Proposed Scheme 0.002 0.62 8.5 0.78 0.12 1.5 Euler-Maruyama 0.005 0.65 9.1 0.75 0.25 2.3 The Proposed Scheme outperforms the Euler-Maruyama method by providing better stability, lower fluctuations, more accurate disease-free equilibrium, and faster computa- tion. It is a more efficient and accurate alternative for solving the stochastic fuzzy SIR-β model. Figure 1: Comparison of proposed and Euler-Maruyama method for stochastic model using d1 = d2 = d3 = d4 = 0.1, µ = 0.1, θ = 0.1, γ = 0.4, c1 = 0.1, Λ = 0.1, β1 = 0.1, α = 0.1, m = 0.1. 5.2 Effect of People’s Resistance on Fuzzy Infective Population Figure 2 illustrates the effect of people’s resistance on the fuzzy infective population by analyzing the behavior of β (transmission rate) over time for different values of α (a resistance-related parameter). The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, γ = 0.1, c1 = 0.1, Λ = 0.1, β1 = 0.1, m = 0.1. The figure compares two fuzzy transmission rates, β1 and β2, under varying values of α (0.1 and 0.3). • Solid Line (β1 at α = 0.1): Represents the transmission rate β1 when the resistance effect is low. At α = 0.1, the transmission rate declines gradually, reducing disease transmission but not drastically. • Dashed Line (β1 at α = 0.3): Represents the case where resistance is increased to α = 0.3. The transmission rate declines more sharply compared to α = 0.1, indicating that higher resistance significantly reduces disease spread. • Dash-Dot Line (β2 at α = 0.1): Shows the behavior of another transmission func- tion β2 at a lower resistance level. The decline is similar to β1 but at a slightly different rate, which could indicate another fuzzy interpretation of disease transmis- sion. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 14 of 24 • Dotted Line (β2 at α = 0.3): Similar to the other case but with higher resistance (α = 0.3), leading to a sharper decline in the transmission rate. Increasing α (resistance factor) decreases β faster, meaning higher resistance leads to lower disease transmission. The difference in β1 and β2 suggests that different interpreta- tions of fuzzy transmission rates exist, but both exhibit the same trend: higher resistance reduces infection spread. This result supports the idea that interventions to increase re- sistance (e.g., vaccination, immunity, or behavioral changes) can effectively control the spread of disease. Figure 2: Effect of people’s resistance on fuzzy infective population: comparison of β1 and β2 under α = 0.1 and α = 0.3. 5.3 Effect of Recovery Rate γ on Fuzzy Infective and Recovered Popula- tions Figure 3 illustrates the recovery rate γ on the dynamics of the infective I and recovered R populations in the fuzzy stochastic SIR-β model. The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, c1 = 0.1, Λ = 0.1, β1 = 0.1, m = 0.1. Left Panel: Shows the infective population I over time for two different recovery rates: γ = 0.1 (slow recovery) and γ = 0.4 (faster recovery). Right Panel: Shows the recovered population R over time for the same recovery rates γ = 0.1 and γ = 0.4. Each panel contains two fuzzy interpretations of I and R, labeled as I1, I2 and R1, R2, with different types of curves. Left Panel (Infective Population I): • Solid and Dashed Lines I1 represent one fuzzy interpretation of the infective popu- lation, while Dash-Dot and Dotted Lines I2 represent another. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 15 of 24 • For γ = 0.1 (slow recovery): The infective population remains higher for a longer time before gradually declining. The peak infection is higher, meaning more indi- viduals are infected due to the slower recovery. • For γ = 0.4 (faster recovery): The infective population declines much faster after reaching its peak. The number of infections is lower, showing that a higher recovery rate effectively reduces disease transmission. Right Panel (Recovered Population R): • Solid and Dashed Lines R1 and Dash-Dot and Dotted Lines R2 represent different fuzzy representations of the recovered population. • For γ = 0.1 (slow recovery): The recovered population increases gradually, reaching a lower peak compared to γ = 0.4. This means that fewer individuals recover quickly, prolonging the epidemic. • For γ = 0.4 (faster recovery): The recovered population increases more rapidly, reaching a higher peak before stabilizing. This indicates that a higher recovery rate leads to more individuals gaining immunity faster, reducing the overall spread of infection. Figure 3: Effect of recovery rate γ on fuzzy infective and recovered populations for γ = 0.1 and γ = 0.4. 5.4 Effect of Incidence Rate Parameter m on Fuzzy Infective Population Figure 4 investigates the effect on the infective population I in the fuzzy stochastic SIR-β model of the incidence rate parameter m. The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, c1 = 0.1, α = 0.1, Λ = 0.1, β1 = 0.1, γ = 0.1. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 16 of 24 The incidence rate parameter m plays a crucial role in controlling disease transmission by incorporating saturation effects in the infection rate. The graph compares two values of m = 0.1 (lower incidence effect) and m = 0.7 (higher incidence effect). Two fuzzy interpretations of the infective population I1 and I2 are represented by different I1 line styles. • Solid and Dashed Lines: represent one fuzzy interpretation of the infective pop- ulation. • Dash-Dot and Dotted Lines I2: represent another fuzzy interpretation of the infective population. • For m = 0.1 (lower incidence effect), the peak of the infective population is higher, suggesting a stronger outbreak. The decline in infection is more gradual, meaning the epidemic lasts longer. • When m = 0.7 (higher incidence effect), the peak infection level is lower, meaning that disease transmission is slowed. The infection decreases more quickly, indicating that the epidemic is under better control. Figure 4: Effect of incidence rate parameter m on fuzzy infective population for m = 0.1 and m = 0.7. 5.5 Comparison of Proposed and NSFD Methods for Susceptible Popu- lation Figure 5 compares the accuracy of the proposed numerical scheme with the existing nonstandard finite difference (NSFD) method in computing the susceptible population S over time and space. The comparison is visualized by plotting the absolute differences M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 17 of 24 between the numerical solutions obtained using the NSFD scheme vs. MATLAB solver pdepe (left plot). The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, c1 = 0.1, α = 0.3, m = 0.1, δ = 0.1, µ = 0.1, Λ = 0.1, γ = 0.1. The comparison is visualized by plotting the absolute differences between the numerical solutions obtained using the NSFD scheme vs. MATLAB solver pdepe (left plot). The right plot shows the absolute error between the proposed scheme and the MATLAB solver pdepe. The error is analyzed over time t and spatial coordinate x. • Left Panel (NSFD vs. MATLAB pdepe Solver): The absolute error |Spdepe − SNSFD| is significantly large, with values reaching close to 0.4. The error is highest at later time steps, indicating that NSFD accumulates numerical inaccuracies over time. The red-colored regions at the boundary suggest that NSFD struggles with maintaining accuracy near boundary conditions. • Right Panel (Proposed Scheme vs. MATLAB pdepe Solver): The absolute error |Spdepe − SProposed| is significantly smaller, remaining below 0.05. The error is distributed more smoothly over time and space, indicating higher numerical stability and accuracy. Unlike NSFD, the proposed scheme avoids large error spikes and remains more consistent throughout the simulation. Figure 5: Comparison of NSFD and proposed scheme with MATLAB pdepe solver for susceptible population S. 5.6 Comparison of Proposed and NSFD Methods for Infected Population Figure 6 compares the accuracy of the proposed numerical scheme with the existing nonstandard finite difference (NSFD) method in computing the infected population I over M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 18 of 24 time and space. The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, c1 = 0.1, α = 0.3, m = 0.1, δ = 0.1, µ = 0.1, Λ = 0.1, γ = 0.1. The comparison is visualized by plotting the absolute differences between the numerical solutions obtained using the NSFD scheme vs. MATLAB solver pdepe (left plot). The right plot shows the absolute error between the proposed scheme and the MATLAB solver pdepe. The absolute error is analyzed over time t and spatial coordinate x. • Left Panel (NSFD vs. MATLAB pdepe Solver): The absolute error |Ipdepe − INSFD| is significantly large, reaching values close to 0.8. The error is highest at later time steps, showing that NSFD accumulates large numerical inaccuracies over time. The red-colored region dominates the plot, indicating that NSFD fails to approxi- mate the infected population accurately. This suggests that the NSFD method lacks numerical stability when modeling the spread of infections. • Right Panel (Proposed Scheme vs. MATLAB pdepe Solver): The absolute error |Ipdepe−IProposed| is significantly smaller, staying below 0.05. The error remains smoothly distributed over space and time, demonstrating higher numerical stability and accuracy. Unlike the NSFD method, the proposed scheme does not accumulate large errors over time, making it more reliable for long-term simulations. Figure 6: Comparison of NSFD and proposed scheme with MATLAB pdepe solver for infected population I. 5.7 Comparison of Proposed and NSFD Methods for Recovered Popula- tion Figure 7 compares the accuracy of the proposed numerical scheme with the existing nonstandard finite difference (NSFD) method in computing the recovered population R M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 19 of 24 over time and space. The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, c1 = 0.1, α = 0.3, m = 0.1, δ = 0.1, µ = 0.1, Λ = 0.1, γ = 0.1. The comparison is visualized by plotting the absolute differences between the numerical solutions obtained using the NSFD scheme vs. MATLAB solver pdepe (left plot). The right plot shows the absolute error between the proposed scheme and the MATLAB solver pdepe. The absolute error is analyzed over time t and spatial coordinate x. • Left Panel (NSFD vs. MATLAB pdepe Solver): The absolute error |Rpdepe − RNSFD| is significantly large, reaching values close to 0.5. The error accumulates over time, indicating that NSFD struggles to track the recovered population dynamics accurately. The red-colored region at the top suggests that NSFD introduces major inaccuracies as time progresses, making long-term simulations unreliable. • Right Panel (Proposed Scheme vs. MATLAB pdepe Solver): The absolute error |Rpdepe−RProposed| is significantly lower, staying below 0.04. The error remains smoothly distributed, showing that the proposed method maintains higher numerical stability. Unlike NSFD, the proposed scheme does not accumulate significant errors over time, making it a better choice for long-term simulations. Figure 7: Comparison of NSFD and proposed scheme with MATLAB pdepe solver for recovered population R. 5.8 Comparison of Proposed and NSFD Methods for β(t, x) (Transmission Rate) Figure 8 compares the accuracy of the proposed numerical scheme with the existing nonstandard finite difference (NSFD) method in computing the transmission rate β(t, x) M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 20 of 24 over time and space. The parameters used in the simulation are: d1 = d2 = d3 = d4 = 0.1, θ = 0.1, c1 = 0.1, α = 0.3, m = 0.1, δ = 0.1, µ = 0.1, Λ = 0.1, γ = 0.1. The comparison is visualized by plotting the absolute differences between the numerical solutions obtained using the NSFD scheme vs. MATLAB solver pdepe (left plot). The right plot shows the absolute error between the proposed scheme and the MATLAB solver pdepe. The absolute error is analyzed over time t and spatial coordinate x. • Left Panel (NSFD vs. MATLAB pdepe Solver): The absolute error |βpdepe − βNSFD| reaches 0.15, indicating a moderate level of error in the NSFD method. The error is not uniform, suggesting that NSFD has difficulty capturing the transmission rate dynamics at certain points in time and space. The error increases over time, showing that NSFD accumulates inaccuracies as the simulation progresses. • Right Panel (Proposed Scheme vs. MATLAB pdepe Solver): The absolute error |βpdepe−βProposed| reaches 0.25, slightly higher than the NSFD error in certain regions. The error distribution is more uniform and smoother, suggesting better numerical stability than NSFD. Despite the slightly higher maximum error, the proposed scheme maintains stability and avoids sharp fluctuations, making it more reliable in long-term simulations. The nonstandard finite difference method used in [51] is also applied here in this study. Since the nonstandard finite difference method is not even first-order accurate in time, the proposed scheme performs better than the nonstandard finite difference scheme. Figure 8: Comparison of NSFD and proposed scheme with MATLAB pdepe solver for transmission rate β(t, x). M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 21 of 24 6. Conclusion With an emphasis on the stochastic fuzzy SIR-β model, including diffusion and a non- linear incidence rate, in this work, we suggested an explicit finite difference scheme for solving fuzzy stochastic partial differential equations. Two-time levels define the developed technique, incorporating a modified time integrator to improve accuracy and stability. To guarantee the dependability of the scheme in solving challenging epidemiological models under uncertainty, we presented a thorough investigation of its stability and consistency in the mean square sense. Applying the suggested approach to the fuzzy stochastic SIR-β model shows how well it captures both stochastic fluctuations and fuzzy uncertainty related to disease spread. We tested, using numerical simulations, the performance of our strategy with an existing non- standard finite difference (NSFD) approach applied to deterministic models. Particularly in cases including stochastic perturbations and fuzzy uncertainty, the findings revealed that our technique beats the NSFD method in terms of accuracy and numerical stability. A numerical scheme has been proposed for solving stochastic time-dependent partial differential equations. A compact scheme is chosen to discretize space-dependent terms. The compact scheme provided sixth-order accuracy on most of the grid points. The scheme has been applied to a mathematical model of SIR-β model. Following are the concluding points: • The proposed scheme performed better than the existing nonstandard finite differ- ence scheme in most cases. • Both groups of infective people declined because of the increase in people’s resistance to change in their social behavior. • Both parts of infected people decline by raising the incidence rate parameter. The results of this work help to progress numerical techniques for uncertain epidemi- ological model solutions. Our method offers a more realistic framework for modeling infectious disease dynamics by combining fuzzy logic with stochastic processes, which can be helpful for public health policy-making. Future studies can expand this work by in- vestigating implicit systems for enhanced stability, higher-order accuracy approaches, and applications to different epidemiological models, including numerous disease compartments and spatial heterogeneity Acknowledgements The authors would like to acknowledge the support of Prince Sultan University for paying the Article Processing Charges (APC) of this publication. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 22 of 24 References [1] W.O. Kermack and A.G. McKendrick. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser., 115:700–721, 1927. [2] A.V. Arundel, E.M. Sterling, J.H. Biggin, and T.D. Sterling. Indirect health effects of relative humidity in indoor environments. Environ. Health Perspect., 65:351–361, 1986. [3] S.M. Minhaz Ud-Dean. Structural explanation for the effect of humidity on persistence of airborne virus: Seasonality of influenza. J. Theoret. Biol., 264:822–829, 2010. [4] M.J. Keeling and P. Rohani. Modeling infectious diseases in human and animals. Princeton University Press, 2008. [5] N. Dexter. Stochastic models of foot and mouth disease in feral pigs in the australian semi-arid rangelands. J. Appl. Ecol., 40:293–306, 2003. [6] M. Liu and K. Wang. Persistence and extinction of a stochastic single-specie model under regime switching in a polluted environment. J. Theor. Biol., 264:934–944, 2010. [7] A. Gray, D. Greenhalgh, X. Mao, and J. Pan. The sis epidemic model with markovian switching. J. Math. Anal. Appl., 394:496–516, 2012. [8] J. Bao and J. Shao. Permanence and extinction of regime-switching predator-prey models. SIAM J. Math. Anal., 48:725–739, 2016. [9] X. Mao. Stochastic differential equations and applications. Elsevier, 2007. [10] N. Bacar and M. Khaladi. On the basic reproduction number in a random environ- ment. J. Math. Biol., 67:1729–1739, 2013. [11] N. Bacar and A. Ed-Darraz. On linear birth-and-death processes in a random envi- ronment. J. Math. Biol., 69:73–90, 2014. [12] S.U. Khan and I. Ali. Application of legendre spectral-collocation method to delay differential and stochastic delay differential equation. AIP Adv., 8:035301, 2008. [13] S.U. Khan, M. Ali, and I. Ali. A spectral collocation method for stochastic volterra integro-differential equations and its error analysis. J. Adv. Differ. Equ., 1:161, 2019. [14] S.U. Khan and I. Ali. Numerical analysis of stochastic sir model by legendre spectral collocation method. Adv. Mech. Eng., 11:1–9, 2019. [15] M. Sohaib. Mathematical modeling and numerical simulation of hiv infection model. Results Appl. Math., 7:100118, 2020. [16] R.D. Attaullah and W. Weera. Galerkin time discretization scheme for the transmis- sion dynamics of hiv infection with non-linear supply rate. J. AIMS Math., 6:11292– 11310, 2022. [17] S. Alyobi and M.F. Yassen. A study on the transmission and dynamical behavior of an hiv/aids epidemic model with a cure rate. AIMS Mathematics, 7(9):17507–17528, 2022. [18] R. Jan and Ş. Yüzbaşı. Dynamical behaviour of hiv infection with the influence of vari- able source term through galerkin method. Chaos, Solitons & Fractals, 152:111429, 2021. [19] N.T. Hieu, N.H. Du, P. Auger, and N.H. Dang. Dynamical behavior of a stochastic sirs epidemic model. Math. Model. Nat. Phenom., 10:56–73, 2015. M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 23 of 24 [20] R. Din and E.A. Algehyne. Mathematical analysis of covid-19 by using sir model with convex incidence rate. Results Phys., 23:103970, 2021. [21] Y.M. Chu, A. Ali, M.A. Khan, S. Islam, and S. Ullah. Dynamics of fractional order covid-19 model with a case study of saudi arabia. Results Phys., 21:103787, 2021. [22] A. Atangana and S.I. Araz. New concept in calculus: Piecewise differential and integral operators. Chaos Solitons Fractals, 145:110638, 2021. [23] A. Atangana and S.I. Araz. Modeling third waves of covid-19 spread with piecewise differential and integral operators. turkey, spain and czechia. medRxiv, 2021. [24] M. Rahman, M. Arfan, K. Shah, and J.F. Gómez-Aguilar. Investigating a non-linear dynamical model of covid-19 disease under fuzzy caputo, random and abc fractional order derivative. Chaos Solitons Fractals, 140:110232, 2020. [25] A. El Koufi, A. Bennar, N. El Koufi, and N. Yousfi. Asymptotic properties of a stochastic siqr epidemic model with lévy jumps and beddington–deangelis incidence rate. Results Phys., 27:104472, 2021. [26] M.B. Matadi. Lie symmetry analysis of stochastic sirs model. Commun. Math. Biol. Neurosci., 23:6–9, 2019. [27] R.R. Wilkinson, F.G. Ball, and K.J. Sharkey. The relationships between message passing, pairwise, kermack–mckendrick and stochastic sir epidemic models. J. Math. Biol., 75:1563–1590, 2017. [28] A. Lahrouz and L. Omari. Extinction and stationary distribution of a stochastic sirs epidemic model with non-linear incidence. Stat. Probab. Lett., 83:960–968, 2013. [29] R. Ullah, M. Farooq, F. Faizullah, M.A. Alghafli, and N. Mlaiki. Fractional stochastic functional differential equations with non-lipschitz condition. AIMS Mathematics, 10(3):7127–7143, 2025. [30] A. Alkhazzan, J. Wang, Y. Nie, S.M.A. Shah, D.K. Almutairi, H. Khan, and J. Alz- abut. Lyapunov-based analysis and worm extinction in wireless networks using stochastic sveir model. Alexandria Engineering Journal, 118:337–353, 2025. [31] 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. [32] V. Capasso. Mathematical Structure of Epidemic Systems, volume 97 of Lecture Notes in Biomathematics. Springer, Berlin, Germany, 1993. [33] B. K. Mishra and A. Prajapati. Spread of malicious objects in computer network: A fuzzy approach. Appl. Appl. Math. Int. J., 8:684–700, 2013. [34] T. W. Liao. A fuzzy multicriteria decision-making method for material selection. J. Manuf. Syst., 15:1–12, 1996. [35] M. S. Arif. A novel explicit scheme for stochastic diffusive sis models with treatment effects. Partial Differential Equations in Applied Mathematics, page 101215, 2025. [36] 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. [37] M. S. Arif, W. Shatanawi, and Y. Nawaz. Stochastic analysis of electro-osmotic flow dynamics in porous media with energy dissipation. International Journal of M. S. Arif, K. Abodayeh, Y. Nawaz / Eur. J. Pure Appl. Math, 18 (3) (2025), 6292 24 of 24 Thermofluids, 27:101172, 2025. [38] M. Saeed, U. Ali, J. Ali, and F. Dayan. Fuzzy soft relative method and its application in decision making problem. Proc. Pak. Acad. Sci. A Phys. Comput. Sci., 57:21–30, 2020. [39] L. D. Barros, M. F. Leite, and R. Bassanezi. The si epidemiological models with a fuzzy transmission parameter. Comput. Math. Appl., 45:1619–1628, 2003. [40] N. R. S. Ortega, P. C. Sallum, and E. Massad. Fuzzy dynamical systems in epidemic modelling. Kybernetes, 29:201–218, 2000. [41] P. K. Mondal, S. Jana, P. Haldar, and T. K. Kar. Dynamical behavior of an epidemic model in a fuzzy transmission. Int. J. Uncertain. Fuzziness Knowl.-Based Syst., 23:651–665, 2015. [42] A. Das and M. Pal. A mathematical study of an imprecise sir epidemic model with treatment control. J. Appl. Math. Comput., 56:477–500, 2018. [43] D. Sadhukhan, L. Sahoo, B. Mondal, and M. Maiti. Food chain model with optimal harvesting in fuzzy environment. J. Appl. Math. Comput., 34:1–18, 2010. [44] C. Li, J. Huang, Y.-H. Chen, and H. Zhao. A fuzzy susceptible-exposed-infected- recovered model based on the confidence index. Int. J. Fuzzy Syst., 23:907–917, 2021. [45] M. Abdy, S. Side, S. Annas, W. Nur, and W. Sanusi. An sir epidemic model for covid- 19 spread with fuzzy parameter: The case of indonesia. Adv. Differ. Equ., pages 1–17, 2021. [46] F. Allehiany et al. Bio-inspired numerical analysis of covid-19 with fuzzy parameters. Comput. Mater. Continua, 72:3213–3229, 2022. [47] R. E. Mickens. Dynamic consistency: A fundamental principle for constructing non- standard finite difference schemes for differential equations. J. Differ. Equ. Appl., 11:645–653, 2005. [48] J. Cresson and F. Pierret. Non standard finite difference scheme preserving dynamical properties. J. Comput. Appl. Math., 303:15–30, 2016. [49] M. Naveed et al. Mathematical analysis of novel coronavirus (2019-ncov) delay pan- demic model. Comput. Mater. Continua, 64:1401–1414, 2020. [50] W. Shatanawi et al. An effective numerical method for the solution of a stochastic coronavirus (2019-ncovid) pandemic model. Comput. Mater. Continua, 66:1121–1137, 2021. [51] S. A. Pasha, Y. Nawaz, and M. S. Arif. On the nonstandard finite difference method for reaction–diffusion models. Chaos, Solitons Fractals, 166:112929, 2023.