EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 5944 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Logistic Growth Epidemiological SEIR Model with Computational and Qualitative Results Kamran1, Sana Maqsood1, Rajermani Thinakaran2, Hasib Khan3,4,∗, Jehad Alzabut3,5 1 Department of Mathematics, Islamia College Peshawar, Peshawar 25120, Khyber Pakhtunkhwa, Pakistan 2 Faculty of Data Science and Information Technology, INTI International University, Malaysia 3 Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833, 11586 Riyadh, Saudi Arabia 4 Department of Mathematics, Shaheed Benazir Bhutto University, Sheringal, Dir Upper, Khyber Pakhtunkhwa, Pakistan 5 Department of Industrial Engineering, OSTIM Technical University, 06374 Ankara, Turkiye Abstract. This work presents a numerical method for solving fractional differential equations arising in mathematical biology. We focus on two models: the fractional logistic growth model, which describes population dynamics, and the fractional SEIR model, which describes the spread of epidemics. We use the fourth-order fractional Runge-Kutta (FRK4) method to approximate the solution of fractional differential equations (FDEs) associated with these models. The Caputo definition of the fractional derivative is used in our models as it is more appropriate for initial value problems involving real-life phenomena. The application of FRK4 provides a stable and accurate numerical solution for both models. The method is validated via numerical experiments, demon- strating efficiency and convergence in controlling fractional-order dynamics, as process innovation. A comparative study with existing numerical techniques highlights the benefits of FRK4 in terms of accuracy and efficiency. 2020 Mathematics Subject Classifications: 92D30, 34A08, 65L06, 37M05 Key Words and Phrases: Fractional Logistic Growth Models, fractional SEIR model, fractional fourth order Runge-Kutta method ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5944 Email addresses: kamran.maths@icp.edu.pk (Kamran), sanamaqsood2002@gmail.com (S. Maqsood), rajermani.thina@newinti.edu.my (T. Thinakaran), hkhan@psu.edu.sa (H. Khan), jalzabut@psu.edu.sa (J. Alzabut), https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 2 of 32 1. Introduction The study of fractional-order problems has now become an area of interest for the research community. Due to the increasing applications of fractional operators, researchers are attracted to the field. Many real-world problems have been modeled and studied using fractional-order operators. Problems involving fractional operators arise in different fields of engineering and other sciences, such as physics, economics, biology, ecology, fluid me- chanics, and electrochemistry [1–4]. Fractional calculation is an efficient mathematical tool to solve different problems in mathematics and engineering. To attract more attention in this area and to validate its applicability and efficiency, in this work, we present a new and recent application of fractional calculus in biological sciences [5, 6]. Recently, fractional calculus has been applied to analyze various nonlinear problems [7, 8]. Since analyti- cal solutions to fractional differential equations are often unavailable, semi-analytical and numerical methods are essential for investigating nonlinear fractional-order problems[9]. Recently, various methods have been developed for solving linear and non-linear dynam- ical systems, such as the Homotopy perturbation method (HPM) [10], Legendre wavelet Tau method[11], the Laplace transform method [12, 13], Adomian decomposition method [14], homotopy analysis transformation method (HATM) [15], the method of operational matrices based on Bernstein polynomials [16], the the extrapolation method [17], and references therein. In literature, numerous numerical methods have been developed for the approximation of linear equations. In [18] the authors have proposed the Adams-Bashforth-Moulton approach for the numerical solutions of nonlinear differential equations of fractional order. The numerical solution of a nonlinear delay differential equation is investigated via the predictor-corrector method in [19]. The authors in [20] have proposed a new predictor- corrector method and compared the results with the fractional Adam’s method for the nonlinear fractional differential equations. Garrapa [21] has discussed the linear stability of the predictor-corrector algorithm for fractional differential equations. Odibat and Momani [22] have developed a new method based on the Fractional Euler’s method and modified the trapezoidal rule by utilizing the generalized Taylor series expansion. Rahimkhani et al. [23] have approximated the solution of a nonlinear pantograph equation of fractional order via the generalized fractional-order Bernoulli wavelet. In [24] the dynamics of fractional order model for COVID-19 have been investigated via an adaptive predictor-corrector and fractional RK4 method. Furthermore, the logistic growth model which is typically the model of population growth is an active area for the research community. The logistic growth model is popular amongst researchers due to its ability to model various biological and social phenomena. The investigation of the solution of the logistic growth model has been considered by many researchers. For example, the authors in [25] have developed a two-stage Runge-Kutta method for fractional logistic growth model. In [26] a spectral tau method is developed for the solution of the fractional logistic growth model. Das et al. [27] discussed the solution of the logistic growth model via the Homotopy perturbation method. Rida et al. [28] have discussed the stability of a fractional SEI model with logistic growth. Other related Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 3 of 32 works on fractional logistic growth model are reported in [29–31] their references. In this paper, we investigate the numerical solutions of fractional logistic growth models and the fractional SEIR model via the fourth-order fractional Runge-Kutta (FRK4) method. This study highlights s a modified SEIR model incorporating logistic growth to better capture real-world epidemic dynamics under constrained population environments. Unlike classical models assuming unlimited growth, the logistic approach consider the saturation effect due to limited resources or interventions. The model studies the transitions between population compartments using a system of nonlinear FDEs. Numerical simulations are conducted to validate theoretical findings and illustrate the impact of parameters such as infection rate and recovery rate on disease transmission. The results show how logistic saturation dampens infection peaks and stabilizes population dynamics over time. This technique suggests a more realistic framework for researchers to evaluate disease control strategies, especially in densely populated areas or during longtime outbreaks. The outline of this paper is as follows: in Section 2, some preliminaries are discussed. In Section 3, the FRK4 method for fractional IVP and fractional SEIR model is illustrated. In section 4, the method is applied to the fractional logistic growth model and the fractional SEIR model. In Section 5, the conclusion is drawn. 2. Preliminaries Definition 1. The Riemann-Liouville (RL) fractional integral of order β > 0 of a function u : (0, ∞) → R is defined by Iβu(t) = 1 Γ(β) ∫ t 0 (t − θ)β−1u(θ)dθ, (2.1) where Γ(·), denotes the gamma function. Definition 2. The RL fractional derivative of order β > 0 of a function u : (0, ∞) → R is defined by Dβu(t) = dm dtm ( 1 Γ(m − β) ∫ t 0 (t − θ)m−β−1u(θ)dθ ) , (2.2) where m = [β] + 1, [·] denotes the integer function. Definition 3. The Caputo fractional derivative of order β > 0 of a function u : (0, ∞) → R is defined by Dβu(t) = 1 Γ(m − β) ∫ t 0 (t − θ)m−β−1 dmu dtm (θ)dθ, (2.3) where m = [β] + 1, [·] denotes the integer function. Here we present some theorems for the existence of a unique solution for the IVP considered in this work. Theorem 1. [32] Assume f ∈ C[R0, R] where R0 = [(t, u : 0 < t < a and |u − u0| ≤ b)] and let |f(t, u)| ≤ M on R0. Then there exists at least one solution to the problem c 0Dβ t u(t) = f (t, u(t)) , u(t0) = u0, 0 < β ≤ 1, 0 ≤ t ≤ η, (2.4) Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 4 of 32 where η = min ( a, [ b M Γ(β + 1)] 1 β ) Theorem 2. [33] Consider the problem (2.4). Let h(µ, u∗(µ)) = f(t − (tβ − µΓ(β + 1)) 1 β , u(t − (tβ − µΓ(β + 1)) 1 β )) and assume that the conditions of Theorem 2.1 hold. Then, a solution u(t) of (2.4) is given by u(t) = u∗( tβ Γ(β + 1)), where u∗(µ) is the solution of the integer order problem du∗(µ) dµ = h(µ, u∗(µ)), with the initial condition u∗(0) = u0 The generalization of the above Theorems 1 and 2 for the system of fractional order differential equations is given by Theorem 3. [33] Let ∥ · ∥ be any norm defined on Rn. Let f ∈ C[R1,Rn], where R1 = [(t, u) : 0 ≤ t ≤ a and∥u − u0∥ ≤ b, f = (f1, f2, ..., fn)T , u = (q1, q2, ..., qn)T , and let ∥f(t, u)∥ ≤ M, R1. Then, there exists at least one solution for the system of fractional differential equations given by c 0Dβ t u(t) = f (t, u(t)) , u(t0) = u0, 0 < β ≤ 1, (2.5) with the initial condition u(0) = u0 on 0 ≤ t ≤ δ, δ = min ( a, [ b M Γ(β + 1)] 1 β ) . 3. Fractional Runge-Kutta Method of Order Four The Fractional Runge-Kutta Method of Order Four (FRK4) is a numerical method designed to solve fractional differential equations (FDEs), extending the classical Runge- Kutta method of order four to handle fractional derivatives, such as those defined in the Riemann-Liouville or Caputo sense. Fractional derivatives, which generalize differentia- tion to non-integer orders, are powerful tools for modeling complex systems with memory effects, such as epidemiological models, viscoelastic materials, and anomalous diffusion processes. The FRK4 method is particularly effective for FDEs because it achieves high accuracy while maintaining numerical stability, making it suitable for solving nonlinear and time-fractional systems. The classical Runge-Kutta method of order four is a widely Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 5 of 32 used numerical technique for solving ordinary differential equations (ODEs) by approxi- mating the solution through a weighted average of four function evaluations per time step. The FRK4 method adapts this approach to FDEs by incorporating the fractional-order derivative, typically denoted by β where 0 < β ≤ 1. Unlike integer-order derivatives, fractional derivatives are non-local, meaning they depend on the history of the solution, which complicates numerical computations. The FRK4 method addresses this by using a generalized Taylor formula tailored to fractional-order systems, as developed by Milici et al. [34]. The FRK4 method offers several advantages: • High Accuracy: It achieves an error of order O(h3β), where h is the step size and β is the fractional order, providing precise solutions for FDEs. • Flexibility: It can handle both Riemann-Liouville and Caputo fractional deriva- tives, which are commonly used in mathematical modeling. • Stability: The method maintains numerical stability for a wide range of FDEs, including nonlinear systems. However, a key limitation is its increased computational cost, particularly when higher accuracy is desired by reducing the step size or increasing the number of evaluations. This cost arises because fractional derivatives require evaluating integrals over the history of the solution, which can be computationally intensive [34]. Despite this, the FRK4 method is often preferred over other numerical methods for FDEs due to its balance of accuracy and robustness. The FRK4 method was initially proposed by Milici et al. [35] to solve FDE and was further refined in [34]. The authors developed a general FRK4 scheme by leveraging the generalized Taylor formula for time-fractional systems. This formula accounts for the non-local nature of fractional derivatives, allowing the method to approximate the solution with a local error function defined as: E(h) = u(tn+1) − u(tn) = ADβu(tn) + BD2βu(tn) + O(h3β), (3.1) where A = hβ Γ(β+1) , B = h2β Γ(2β+1) , Dβ denotes the fractional derivative of order β, O(·) is the Bachmann-Landau notation for higher-order terms [34], and Γ(·) is the Gamma function, which generalizes factorials to real and complex numbers. This error function quantifies the truncation error per step, showing that the accuracy of the method improves with smaller step sizes but at the cost of more computations. 3.1. Differential Equation of Fractional Order To illustrate the FRK4 method, consider an initial value problem (IVP) for a fractional- order differential equation of the form: C 0 Dβ t u(t) = f(t, u(t)), u(t0) = u0, u(t) ∈ Cp+1([t0, t0 + T ]), 0 < β ≤ 1, (3.2) Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 6 of 32 where C 0 Dβ t is the Caputo fractional derivative of order β, defined as: C 0 Dβ t u(t) = 1 Γ(1 − β) ∫ t 0 (t − τ)−βu′(τ)dτ, f(t, u(t)) is a nonlinear function governing the dynamics, u0 is the initial condition, Cp+1 denotes the space of functions with p + 1-th order continuous derivatives, and T is the time interval. The Caputo derivative is preferred in many applications because it allows initial conditions to be specified in terms of integer-order derivatives, which aligns with physical interpretations (e.g., initial population in epidemiological models). 3.2. FRK4 Numerical Scheme The FRK4 method approximates the solution of (3.2) at discrete time points tn = t0 + nh, where h is the size of the step. The numerical solution at the next time step, un+1, is computed as: un+1 = un + A 6 (K1 + 2K2 + 2K3 + K4), (3.3) where A = hβ Γ(β+1) and the stage values K1, K2, K3, K4 are given by: K1 = f(tn, un), K2 = f ( tn + A 2 , un + A 2 K1 ) , K3 = f ( tn + A 2 , un + A 2 K2 ) , K4 = f(tn + A, un + AK3). These stage values approximate the function f(t, u(t)) at intermediate points, weighted to account for the non-local behavior of the fractional derivative. The coefficients (1, 2, 2, 1) mirror the classical Runge-Kutta method, but the time increments and arguments are scaled by A, which incorporates the fractional order β. To apply the FRK4 method, the following steps are required: The Basic steps of the proposed numerical method are presented in Algorithm 1. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 7 of 32 Algorithm 1 : Algorithm for the FRK4 method 1: Input: initial condition u0, initial time t0, final time T , step size h, function f(t, u), fractional order β. 2: Outputs: arrays t (time points) and u (solutions). 3: Step (a): Discretize time compute n = T −t0 h and set tn = t0 + nh. 4: Step (b): Initialize u[0] = u0 . 5: Step (c): Iterate loop over n, compute A, evaluate stages K1 to K4 and update un+1. 6: Step (d): Return output t and u. 3.3. Limitations and Challenges While the FRK4 method is robust, its computational cost is a significant drawback. Evaluation of K1 to K4 requires multiple function calls per step, and the non-local nature of the fractional derivative increases memory requirements, especially for small h or large T . Additionally, the method’s accuracy depends on the smoothness of f(t, u(t)), and discontinuities or singularities in the solution can reduce performance. 3.4. Fractional SEIR Model The SEIR model is a cornerstone of mathematical epidemiology, used to describe the spread of infectious diseases by dividing a population into four compartments: Susceptible (S), Exposed (E), Infected (I ), and Recovered (R). In this section, we employ the Frac- tional Runge-Kutta Method of Order Four (FRK4) to numerically solve a fractional-order SEIR model, which incorporates Caputo fractional derivatives of order 0 < β < 1 to cap- ture memory effects and non-local dynamics in disease transmission. Fractional derivatives are particularly suited for epidemiological models because they account for historical de- pendencies, such as prolonged incubation periods or persistent immunity, which standard integer-order derivatives cannot model effectively. The fractional SEIR model is defined as: c 0Dβ t S = bN − pbE − qbI − r SI N − d(N)S, c 0Dβ t E = pbE + qbI + r SI N − d(N)E, c 0Dβ t I = ηE − θI − γI − d(N)I, c 0Dβ t R = γI − d(N)R. (3.4) S(0) = S0, E(0) = E0, I(0) = I0, R(0) = R0, 0 < β < 1. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 8 of 32 Here, N = S + E + I + R is the total population, assumed constant for simplicity. The parameters have the following epidemiological interpretations: • b: Birth rate (new individuals entering the susceptible population). • p: Probability of exposure from contact with exposed individuals. • q: Probability of exposure from contact with infected individuals. • r: Transmission rate due to interactions between susceptible and infected individuals. • d(N): Natural death rate, potentially density-dependent. • η: Rate at which exposed individuals become infectious. • θ: Disease-induced mortality rate. • γ: Recovery rate of infected individuals. The fractional derivative C 0 Dβ t , introduces memory effects, allowing the model to reflect phenomena like delayed disease progression or long-term immunity, with β = 1 recovering the classical SEIR model. To apply the FRK4 method, we rewrite the system by defining the right-hand side functions: f1(t, S, E, I, R) = bN − pbE − qbI − r SI N − d(N)S, f2(t, S, E, I, R) = pbE + qbI + r SI N − d(N)E, f3(t, S, E, I, R) = ηE − θI − γI − d(N)I, f4(t, S, E, I, R) = γI − d(N)R. Thus, the SEIR model becomes: c 0Dβ t S = f1(t, S, E, I, R) c 0Dβ t E = f2(t, S, E, I, R) c 0Dβ t I = f3(t, S, E, I, R) c 0Dβ t R = f4(t, S, E, I, R) Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 9 of 32 The FRK4 method approximates the solution at discrete time points tn = t0 + nh, where h is the step size, by computing: Sn+1 = Sn + A 6 (K1S + 2K2S + 2K3S + K4S), En+1 = En + A 6 (K1E + 2K2E + 2K3E + K4E), In+1 = In + A 6 (K1I + 2K2I + 2K3I + K4I), Rn+1 = Rn + A 6 (K1R + 2K2R + 2K3R + K4R), where A = hβ Γ(β+1) , and the stage values are: K1S = f1(tn, Sn, En, In, Rn), K1E = f2(tn, Sn, En, In, Rn), K1I = f3(tn, Sn, En, In, Rn), K1R = f4(tn, Sn, En, In, Rn), K2S = f1 ( tn + A 2 , Sn + A 2 K1S , En + A 2 K1E , In + A 2 K1I , Rn + A 2 K1R ) , K2E = f2 ( tn + A 2 , Sn + A 2 K1S , En + A 2 K1E , In + A 2 K1I , Rn + A 2 K1R ) , K2I = f3 ( tn + A 2 , Sn + A 2 K1S , En + A 2 K1E , In + A 2 K1I , Rn + A 2 K1R ) , K2R = f4 ( tn + A 2 , Sn + A 2 K1S , En + A 2 K1E , In + A 2 K1I , Rn + A 2 K1R ) , K3S = f1 ( tn + A 2 , Sn + A 2 K2S , En + A 2 K2E , In + A 2 K2I , Rn + A 2 K2R ) , K3E = f2 ( tn + A 2 , Sn + A 2 K2S , En + A 2 K2E , In + A 2 K2I , Rn + A 2 K2R ) , K3I = f3 ( tn + A 2 , Sn + A 2 K2S , En + A 2 K2E , In + A 2 K2I , Rn + A 2 K2R ) , K3R = f4 ( tn + A 2 , Sn + A 2 K2S , En + A 2 K2E , In + A 2 K2I , Rn + A 2 K2R ) , K4S = f1(tn + A, Sn + AK3S , En + AK3E , In + AK3I , Rn + AK3R), K4E = f2(tn + A, Sn + AK3S , En + AK3E , In + AK3I , Rn + AK3R), K4I = f3(tn + A, Sn + AK3S , En + AK3E, In + AK3I , Rn + AK3R), K4R = f4(tn + A, Sn + AK3S , En + AK3E , In + AK3I , Rn + AK3R). The FRK4 method is applied to each compartment simultaneously, updating S, E, I, and R at each time step. Further details on numerical methods with applications to the Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 10 of 32 FDEs can be found in [34, 36–41]. 4. Results and Discussion This section is devoted to the numerical investigation of different fractional order models by employing the FRK4 method. 4.0.1. Problem 1 Here we consider the logistic growth model of fractional order of the form [25] c 0Dβ t u(t) = γu(t) ( 1 − u(t) M ) (4.1) u(t0) = u0, where u0 is the initial density , and γ is the intrinsic growth rate of the population, and M is the carrying capacity. The exact solution of the problem for β = 1 is given by u(t) = Mu0 u0+(M−u0)e−γt . Using Theorem 2 the exact solution of the given problem for noninteger values becomes u(t) = Mu0 u0 + (M − u0)e −γ ( tβ Γ(β+1) ) . The numerical results obtained using the Fractional Runge-Kutta Method of Order Four (FRK4) for solving a fractional differential equation are presented in Tables 1 to 4, with corresponding plots in Figures 1a to 6b. The parameters include the initial condition u0, fractional order β, step size h, growth parameter γ, and a model-specific parameter M . Tables 1 and 2 show results for M = 10, γ = 0.5, u0 = 20, with h = 0.01 and h = 0.005, respectively, for β = 0.96 and β = 1. Tables 3 and 4 present results for M = 1, γ = 0.5, with u0 = 20, h = 0.01 (Table 3) and u0 = 10, h = 0.005 (Table 4), again for β = 0.96 and β = 1. The results demonstrate that the FRK4 method achieves high accuracy, with errors decreasing as the step size h is reduced. Further, these plots demonstrate the sensitivity of the model to β. As we see that for a lower value of β (e.g., 0.96), the population approaches the carrying capacity more slowly, reflecting stronger memory effects that delay growth or decline. Figures 1a and 1b illustrate the exact and approximate solutions for M = 10, u0 = 20, γ = 0.5, h = 0.01, with β = 1 and β = 0.96, respectively, showing close agreement between numerical and analytical solutions. Similarly, Figures 2a and 2b depict solutions for M = 1, u0 = 20, γ = 0.5, h = 0.005, for β = 1 and β = 0.96, while Figures 4a and 4b show the same for u0 = 10. Figures 3a and 3b appear to duplicate the conditions of Figures 2a and 2b, confirming consistent results. Figures 5a and 5b plot numerical solutions for varying β with M = 1, γ = 0.5, h = 0.005, and u0 = 10 and u0 = 20, respectively, highlighting the impact of the fractional order on solution behavior. Absolute error plots for problem 1, with M = 10, γ = 0.5, h = 0.005, β = 1, Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 11 of 32 and u0 = 20 (Figure 6a) and u0 = 10 (Figure 6b), demonstrate low errors underscoring the method’s precision. The FRK4 method’s efficiency is evident across all tables and figures, consistently outperforming the method proposed in [25] in terms of accuracy. The FRK4 method produced significantly more accurate results, particularly for smaller step sizes and fractional orders β < 1, due to its ability to handle the non-local nature of fractional derivatives. The observed trend of improved accuracy with decreasing h aligns with the method’s theoretical error bound, making it a robust choice for solving fractional differential equations for further details on the FRK4 method, see [34]. Table 1: The numerical results for problem 1 using FRK4 method with M = 10, γ = 0.5, u0 = 20, h = 0.01. 1* β = 0.96 β = 1 t uapprox uexact errorabs [25] uapprox uexact errorabs [25] 0 20.0000 20.0000 0 0 20.0000 20.0000 0 0 0.1 18.8778 19.0699 9.0833×10−2 0.1553 19.0678 19.0699 1.2496×10−10 2.3897×10−5 0.2 17.9316 18.2621 2.0004×10−1 0.2680 18.2586 18.2621 1.9640×10−10 3.9362×10−5 0.3 17.1245 17.5548 2.9312×10−1 0.3500 17.5502 17.5548 2.3539×10−10 4.9190×10−5 0.4 16.4290 16.9309 3.6881×10−1 0.4093 16.9257 16.9309 2.5442×10−10 5.5198×10−5 0.5 15.8246 16.3773 4.2913×10−1 0.4518 16.3717 16.3773 2.6107×10−10 5.8593×10−5 0.6 15.2954 15.8833 4.7654×10−1 0.4816 15.8774 15.8833 2.6005×10−10 6.0190×10−5 0.7 14.8290 15.4403 5.1326×10−1 0.5017 15.4342 15.4403 2.5433×10−10 6.0546×10−5 0.8 14.4156 15.0412 5.4120×10−1 0.5144 15.0351 15.0412 2.4581×10−10 6.0048×10−5 0.9 14.0472 14.6803 5.6194×10−1 0.5213 14.6741 14.6803 2.3571×10−10 5.8968×10−5 1 13.7176 14.3527 5.7676×10−1 0.5238 14.3466 14.3527 2.2482×10−10 5.7498×10−5 Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 12 of 32 Table 2: The numerical results for problem 1 using FRK4 method with h = 0.005, M = 10, γ = 0.5, u0 = 20. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 20.0000 20.0000 0 20.0000 20.0000 0 0.1 18.8505 19.0699 1.1957×10−1 19.0689 19.0699 7.7982×10−12 0.2 17.8853 18.2621 2.4873×10−1 18.2604 18.2621 1.2253×10−11 0.3 17.0649 17.5548 3.5568×10−1 17.5525 17.5548 1.4683×10−11 0.4 16.3602 16.9309 4.4093×10−1 16.9283 16.9309 1.5863×10−11 0.5 15.7496 16.3773 5.0771×10−1 16.3745 16.3773 1.6275×10−11 0.6 15.2163 15.8833 5.5930×10−1 15.8804 15.8833 1.6216×10−11 0.7 14.7474 15.4403 5.9854×10−1 15.4373 15.4403 1.5858×10−11 0.8 14.3327 15.0412 6.2776×10−1 15.0381 15.0412 1.5326×10−11 0.9 13.9639 14.6803 6.4883×10−1 14.6772 14.6803 1.4703×10−11 1 13.6345 14.3527 6.6328×10−1 14.3496 14.3527 1.4026×10−11 Table 3: The numerical results for problem 1 using FRK4 method with M = 1, γ = 0.5, h = 0.01, u0 = 20. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 20.0000 20.0000 0 20.0000 20.0000 0 0.1 9.2732 10.3808 4.4651×100 10.2652 10.3808 5.3672×10−6 0.2 6.1926 7.1223 4.9691×10−1 7.0410 7.1223 2.6602×10−6 0.3 4.7241 5.4846 4.6844×10−1 5.4270 5.4846 1.5369×10−6 0.4 3.8645 4.5003 4.2969×10−1 4.4576 4.5003 9.9278×10−7 0.5 3.3006 3.8441 3.9314×10−1 3.8111 3.8441 6.9170×10−7 0.6 2.9027 3.3758 3.6096×10−1 3.3495 3.3758 5.0846×10−7 0.7 2.6073 3.0253 3.3304×10−1 3.0038 3.0253 3.8888×10−7 0.8 2.3797 2.7533 3.0880×10−1 2.7355 2.7533 3.0662×10−7 0.9 2.1992 2.5364 2.8762×10−1 2.5213 2.5364 2.4764×10−7 1 2.0528 2.3596 2.6897×10−1 2.3467 2.3596 2.0394×10−7 Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 13 of 32 Table 4: The numerical results for problem 1 using FRK4 method with M = 1, u0 = 10, γ = 0.5, h = 0.005. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 10.0000 10.0000 0 10.0000 10.0000 0 0.1 6.4423 6.9496 2.6127×10−1 6.9360 6.9496 1.3360×10−8 0.2 4.8325 5.3866 3.4761×10−1 5.3733 5.3866 9.6239×10−9 0.3 3.9152 4.4373 3.6279×10−1 4.4259 4.4373 6.6909×10−9 0.4 3.3234 3.8002 3.5485×10−1 3.7907 3.8002 4.8216×10−9 0.5 2.9105 3.3436 3.3838×10−1 3.3356 3.3436 3.6101×10−9 0.6 2.6066 3.0006 3.2152×10−1 2.9938 3.0006 2.7923×10−9 0.7 2.3739 2.7339 3.0383×10−1 2.7280 2.7339 2.2180×10−9 0.8 2.1903 2.5207 2.8700×10−1 2.5157 2.5207 1.8008×10−9 0.9 2.0421 2.3467 2.7129×10−1 2.3422 2.3467 1.4887×10−9 1 1.9200 2.2020 2.5676×10−1 2.1981 2.2020 1.2494×10−9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 14 15 16 17 18 19 20 u (t ) Frac-RK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 13 14 15 16 17 18 19 20 u (t ) FRK4 Exact (b) Figure 1: (a) The approximate and exact solution of problem 1 using FRK4 method with M = 10, u0 = 20, γ = 0.5, h = 0.01, β = 1. We see that the approximate solution closely follows the exact solution. (b) The approximate and exact solution of problem 1 using FRK4 method with M = 10, u0 = 20, γ = 0.5, h = 0.01, β = 0.96. The approximate solution for β = 0.96 deviates slightly, reflecting the memory effects introduced by the fractional derivative. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 14 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 14 15 16 17 18 19 20 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 13 14 15 16 17 18 19 20 u (t ) FRK4 Exact (b) Figure 2: (a) The approximate and exact solution of problem 1 using FRK4 method with M = 10, u0 = 20, γ = 0.5, h = 0.005, β = 1. Here we observed that the smaller step size h = 0.005 improves the accuracy. (b) The approximate and exact solutions of problem 1 using FRK4 method with M = 10, u0 = 20, γ = 0.5, h = 0.005, β = 0.96. Here, we observe that the numerical solution deviates, but the results are still acceptable. The figures highlight that reducing h enhances the method’s ability to approximate the exact solution. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 2 4 6 8 10 12 14 16 18 20 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 2 4 6 8 10 12 14 16 18 20 u (t ) Frac-RK4 Exact (b) Figure 3: (a) The approximate and exact solution of problem 1 using FRK4 method with M = 1, u0 = 20, γ = 0.5, h = 0.005, β = 1. (b) The approximate and exact solution of problem 1 using FRK4 method with M = 1, u0 = 20, γ = 0.5, h = 0.005, β = 0.96. The lower carrying capacity M = 1 results in a faster decline from the initial condition u0 = 20. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 15 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 2 3 4 5 6 7 8 9 10 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 1 2 3 4 5 6 7 8 9 10 u (t ) FRK4 Exact (b) Figure 4: (a) The approximate and exact solution of problem 1 using FRK4 method with M = 1, u0 = 10, γ = 0.5, h = 0.005, β = 1. (b) The approximate and exact solutions of problem 1 using FRK4 method with M = 1, u0 = 10, γ = 0.5, h = 0.005, β = 0.96. The initial condition u0 = 10 is closer to the carrying capacity. The figures show that the approximate solutions track the exact solutions well, but the slower dynamics for beta = 0.96 suggest that memory effects are significant even with a lower initial population, impacting predictions of population stabilization. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 1 2 3 4 5 6 7 e r r o r a b s 10 -14 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 e r r o r a b s 10 -14 (b) Figure 5: (a) The errorabs of FRK4 for problem 1 with M = 1, u0 = 10, γ = 0.5, h = 0.0001, for β = 1. (b) The errorabs of FRK4 for problem 1 with M = 10, u0 = 20, γ = 0.5, h = 0.0001, for β = 1. Both figures confirm the FRK4 method’s high accuracy for β = 1. The smaller errors for M = 10 compared to M = 1 suggest that the method performs better for higher carrying capacities, possibly due to less rapid changes in the solution. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 16 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 e r r o r a b s 10 -11 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 e r r o r a b s 10 -8 (b) Figure 6: (a) The errorabs of FRK4 for problem 1 with M = 10, u0 = 20, γ = 0.5, h = 0.005 β = 1. (b) The errorabs of FRK4 for problem 1 with M = 1, u0 = 10, γ = 0.5, h = 0.005 β = 1.Here we see that the errors are comparatively larger for a large step. 4.0.2. Problem 2 Here we consider exponential growth of Malthus model of the form [5] c 0Dβ t u(t) = ku(t), (4.2) with the initial condition u0 = 1. The exact solution of the problem is u(t) = u0Eβ(ktβ). The results obtained using the FRK4 method for the fractional orders β = 1 and β = 0.96, with parameters u0 = 1, h = 0.01, and k = 1, are presented in Table 5. Figures 7a and 7b illustrate the numerical and exact solutions for β = 1 and β = 0.96, respectively, under the same parameter settings. Figure 8a shows the approximate solutions for different values of β with u0 = 1, k = 1, and step size h = 0.1. Similarly, Figure 8b presents the approximate solutions for various β values with a finer step size h = 0.01. The absolute error errorabs for Problem 2 with β = 1 is plotted in Figure 9a for h = 0.0001, and in Figure 9b for h = 0.0005. Overall, the numerical results show good agreement with those obtained using the method proposed in [5], thereby validating the accuracy and reliability of the FRK4 method for solving fractional differential equations. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 17 of 32 Table 5: The numerical results for problem 2 using FRK4 method with u0 = 1, h = 0.01, k = 1. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 1.0000 1.0000 0 1.0000 1.0000 0 0.1 1.1293 1.1182 1.1416×10−2 1.1050 1.1052 3.6526×10−13 0.2 1.2734 1.2433 3.0692×10−2 1.2210 1.2214 7.6850×10−13 0.3 1.4334 1.3799 5.4445×10−2 1.3493 1.3499 1.2146×10−12 0.4 1.6109 1.5299 8.2293×10−2 1.4910 1.4918 1.7075×10−12 0.5 1.8076 1.6950 1.1426×10−1 1.6476 1.6487 2.2524×10−12 0.6 2.0255 1.8771 1.5056×10−1 1.8208 1.8221 2.8542×10−12 0.7 2.2667 2.0779 1.9147×10−1 2.0121 2.0138 3.5194×10−12 0.8 2.5338 2.2996 2.3737×10−1 2.2235 2.2255 4.2548×10−12 0.9 2.8292 2.5444 2.8871×10−1 2.4572 2.4596 5.0684×10−12 1 3.1561 2.8146 3.4599×10−1 2.7154 2.7183 5.9655×10−12 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 1 1.5 2 2.5 3 3.5 u (t ) FRK4 Exact (b) Figure 7: (a) The approximate and exact solution of problem 2 using FRK4 method with u0 = 1, k = 1, h = 0.01, β = 1. We see that the approximate solution is in good agreement with the exact solution. (b) The approximate and exact solution of problem 2 using the FRK4 method with u0 = 1, k = 1, h = 0.01, β = 0.96. The figure shows a slight deviation from the exact solution, suggesting that lower β values introduce memory effects that temper exponential growth, relevant for modeling populations with historical constraints. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 18 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 2 4 6 8 10 12 14 u (t ) =0.5 =0.6 =0.7 =0.8 =0.9 =1 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 5 10 15 20 25 30 35 40 u (t ) =0.5 =0.6 =0.7 =0.8 =0.9 =1 (b) Figure 8: (a) The approximate solutions of problem 2 using FRK4 method with u0 = 1, k = 1, h = 0.1 for different fractional orders β. Note that the lower the order of the derivative, the greater the rate of variation. (b) The approximate solutions of problem 2 using FRK4 method with u0 = 1, k = 1, h = 0.01 for different fractional orders β. The smaller step size in this figure results in smoother curves, indicating improved accuracy. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 e r r o r a b s 10 -14 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 e r r o r a b s 10 -15 (b) Figure 9: (a) The errorabs of FRK4 for problem 2 with u0 = 1, k = 1, h = 0.0001 for β = 1. (b) The errorabs of FRK4 for problem 2 with u0 = 1, k = 1, h = 0.0005 for β = 1. The smaller h in Figure (a) yields slightly lower errors; however, the stable, low errors indicate that FRK4 is highly effective for exponential growth models, particularly when high accuracy is needed for long-term predictions. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 19 of 32 4.1. Problem 3 Here we consider the quadratic logistic growth model of the form c 0Dβ t u(t) = 1 2u(t)(1 − u(t)), (4.3) with the initial condition u0 = 1 4 . The exact solution of the problem for β = 1 is given as u(t) = e t 2 e t 2 +3 . Using Theorem 2the the exact solution of the given problem for non integer value of β will be of the form u(t) = e ( tβ 2Γ(β+1) ) e ( tβ 2Γ(β+1) ) + 3 The results obtained using the FRK4 method for fractional orders β = 1 and β = 0.96, with u0 = 0.25 and h = 0.1, are presented in Table 6. Additionally, Table 7 presents the results for the same fractional orders using a finer step size h = 0.01, with u0 = 0.25. Figures 10a and 10b show the plots of the numerical and exact solutions for β = 1 and β = 0.96, respectively, using u0 = 0.25 and h = 0.1. Similarly, Figures 11a and 11b depict the corresponding plots for h = 0.01. Figure 12a displays the approximate solutions for various values of β with u0 = 0.25 and h = 0.1, while Figure 12b shows the same for h = 0.01. The absolute error errorabs for Problem 3 with β = 1 is illustrated in Figure 13a for h = 0.01, and in Figure 13b for a finer step size h = 0.0001. Table 6: The numerical results for problem 3 using FRK4 method with u0 = 0.25, h = 0.1. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 0.2500 0.2500 0 0.2500 0.2500 0 0.1 0.2605 0.2595 1.9586×10−10 0.2595 0.2595 1.1411×10−10 0.2 0.2714 0.2692 6.0238×10−4 0.2691 0.2692 2.2664×10−10 0.3 0.2825 0.2792 1.4551×10−3 0.2790 0.2792 3.3717×10−10 0.4 0.2939 0.2893 2.4901×10−3 0.2892 0.2893 4.4531×10−10 0.5 0.3055 0.2997 3.6755×10−3 0.2995 0.2997 5.5068×10−10 0.6 0.3174 0.3103 4.9920×10−3 0.3101 0.3103 6.5295×10−10 0.7 0.3296 0.3211 6.4259×10−3 0.3208 0.3211 7.5182×10−10 0.8 0.3420 0.3321 7.9661×10−3 0.3318 0.3321 8.4701×10−10 0.9 0.3546 0.3433 9.6030×10−3 0.3430 0.3433 9.3832×10−10 1 0.3674 0.3547 1.1327×10−2 0.3543 0.3547 1.0256×10−9 Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 20 of 32 Table 7: The numerical results for problem 3 using FRK4 method with u0 = 0.25, h = 0.01. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 0.2500 0.2500 0 0.2500 0.2500 0 0.1 0.2616 0.2595 1.0373×10−3 0.2595 0.2595 1.1546×10−14 0.2 0.2736 0.2692 2.7344×10−3 0.2692 0.2692 2.2871×10−14 0.3 0.2859 0.2792 4.7370×10−3 0.2792 0.2792 3.4084×10−14 0.4 0.2985 0.2893 6.9745×10−3 0.2893 0.2893 4.4964×10−14 0.5 0.3115 0.2997 9.4118×10−3 0.2997 0.2997 5.5567×10−14 0.6 0.3247 0.3103 1.2026×10−2 0.3103 0.3103 6.5892×10−14 0.7 0.3383 0.3211 1.4799×10−2 0.3211 0.3211 7.5884×10−14 0.8 0.3521 0.3321 1.7715×10−2 0.3321 0.3321 8.5487×10−14 0.9 0.3661 0.3433 2.0760×10−2 0.3433 0.3433 9.4702×10−14 1 0.3804 0.3547 2.3917×10−2 0.3546 0.3547 1.0347×10−13 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.24 0.26 0.28 0.3 0.32 0.34 0.36 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 u (t ) FRK4 Exact (b) Figure 10: (a) The exact and numerical solution of problem 3 using FRK4 method with u0 = 0.25, h = 0.1, β = 1. The numerical solution shows good agreement with the exact solution. (b) The exact and numerical solution of problem 3 using the FRK4 method with u0 = 0.25, h = 0.1, β = 0.96. The slight lag observed in the numerical solution suggests that decreasing β delays convergence to the equilibrium point, which is characteristic of populations under limited growth conditions. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 21 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.24 0.26 0.28 0.3 0.32 0.34 0.36 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.25 0.3 0.35 0.4 u (t ) FRK4 Exact (b) Figure 11: (a) The approximate and exact solution of problem 3 using the FRK4 method with u0 = 0.25, h = 0.01, β = 1. The numerical results are in excellent agreement with the exact analytical solution, demonstrating the accuracy of the proposed method. (b) The approximate and exact solutions of problem 3 using the FRK4 method with u0 = 0.25, h = 0.01, β = 0.96. With a smaller step size of h = 0.01, the approximate solution shows slightly greater deviation. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 u (t ) =0.5 =0.6 =0.7 =0.8 =0.9 =1 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u (t ) =0.5 =0.6 =0.7 =0.8 =0.9 =1 (b) Figure 12: (a) The numerical solutions of problem 3 using FRK4 method with u0 = 0.25, h = 0.1 using different fractional orders β. (b) The numerical solutions of problem 3 using FRK4 method with u0 = 0.25, h = 0.01 using different fractional orders β. Lower β values result in slower growth toward the equilibrium. The smaller h in Figure (b) produces smoother solutions, indicating higher accuracy. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 22 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.2 0.4 0.6 0.8 1 1.2 e r r o r a b s 10 -13 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 e r r o r a b s 10 -15 (b) Figure 13: (a) The errorabs of FRK4 for problem 3 with u0 = 0.25, h = 0.01 and orders β = 1. (b) The errorabs of FRK4 for problem 3 with u0 = 0.25, h = 0.0001 and orders β = 1. The profiles indicate that the error further reduced with h = 0.0001, demonstrating the robustness of the method. 4.2. Problem 4 Here we consider the quadratic fractional logistic growth model. c 0Dβ t u(t) = 1 2u(t)(1 − u(t)), (4.4) with the initial condition u0 = 0.5. The exact solution of the problem for β = 1 is given as u(t) = e t 2 e t 2 +1 . Using Theorem 4.5 the exact solution of the given problem for non-integer value of β will be of the form u(t) = e ( tβ 2Γ(β+1) ) e ( tβ 2Γ(β+1) ) + 1 The results obtained using the FRK4 method for fractional orders β = 1 and β = 0.96, with u0 = 0.5 and h = 0.1, are presented in Table 8. Additionally, Table 9 presents the results for the same fractional orders with u0 = 0.5 and a finer step size h = 0.01. Figures 14a and 14b show the plots of numerical and exact solutions for β = 1 and β = 0.96, respectively, using u0 = 0.5 and h = 0.1. Similarly, Figures 15a and 15b depict the corresponding plots for h = 0.01. Figure 16a illustrates the approximate solutions for various values of β with u0 = 0.5 and h = 0.01, while Figure 16b shows the same for h = 0.1. The absolute error errorabs for Problem 4 with β = 1 is plotted in Figure 17a Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 23 of 32 for h = 0.01, and for Problem 3 in Figure 17b with h = 0.0001. Table 8: The numerical results for problem 4 using FRK4 method with u0 = 0.5, h = 0.1. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 0.5000 0.5000 0 0.5000 0.5000 0 0.1 0.5139 0.5125 7.0019×10−11 0.5125 0.5125 4.0700×10−11 0.2 0.5278 0.5250 7.5967×10−4 0.5250 0.5250 8.1672×10−11 0.3 0.5417 0.5374 1.7850×10−3 0.5374 0.5374 1.2327×10−10 0.4 0.5555 0.5498 2.9711×10−3 0.5499 0.5498 1.6584×10−10 0.5 0.5693 0.5622 4.2652×10−3 0.5622 0.5622 2.0971×10−10 0.6 0.5829 0.5744 5.6335×10−3 0.5745 0.5744 2.5520×10−10 0.7 0.5964 0.5866 7.0516×10−3 0.5867 0.5866 3.0260×10−10 0.8 0.6098 0.5987 8.5003×10−3 0.5988 0.5987 3.5216×10−10 0.9 0.6230 0.6106 9.9637×10−3 0.6107 0.6106 4.0412×10−10 1 0.6360 0.6225 1.1428×10−2 0.6226 0.6225 4.5866×10−10 Table 9: The numerical results for problem 4 using FRK4 method with u0 = 0.5, h = 0.01. 1* β = 0.96 β = 1 t uapprox uexact errorabs uapprox uexact errorabs 0 0.5000 0.5000 0 0.5000 0.5000 0 0.1 0.5153 0.5125 1.3430×10−3 0.5125 0.5125 3.8858×10−15 0.2 0.5305 0.5250 3.4389×10−3 0.5250 0.5250 7.8826×10−15 0.3 0.5457 0.5374 5.7867×10−3 0.5374 0.5374 1.1879×10−14 0.4 0.5608 0.5498 8.2751×10−3 0.5498 0.5498 1.6098×10−14 0.5 0.5758 0.5622 1.0845×10−2 0.5622 0.5622 2.0428×10−14 0.6 0.5907 0.5744 1.3456×10−2 0.5744 0.5744 2.5091×10−14 0.7 0.6053 0.5866 1.6078×10−2 0.5866 0.5866 2.9754×10−14 0.8 0.6199 0.5987 1.8687×10−2 0.5987 0.5987 3.4750×10−14 0.9 0.6341 0.6106 2.1262×10−2 0.6106 0.6106 3.9746×10−14 1 0.6482 0.6225 2.3785×10−2 0.6225 0.6225 4.5408×10−14 Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 24 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 u (t ) FRK4 Exact (b) Figure 14: (a) The approximate and exact solutions of problem 4 using the FRK4 method with u0 = 0.5, h = 0.1, β = 1. The approximate solution is nearly indistinguishable from the exact solution. (b) The approximate and exact solutions of problem 4 using the FRK4 method with u0 = 0.5, h = 0.1, β = 0.96. The Figure shows a minor deviation for beta = 0.96, reflecting memory effects that slow growth, highlighting the FRK4’s effectiveness in capturing fractional dynamics. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 u (t ) FRK4 Exact (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.5 0.55 0.6 0.65 u (t ) FRK4 Exact (b) Figure 15: (a) The approximate and exact solutions of problem 4 using the FRK4 method with u0 = 0.5, h = 0.01, β = 1. The approximate and exact solutions are in good agreement. (b) The approximate and exact solutions of problem 4 using the FRK4 method with u0 = 0.5, h = 0.01, β = 0.96. A smaller step size h leads to a slightly greater deviation in the approximate solution, as observed in the numerical results. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 25 of 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 u (t ) =0.5 =0.6 =0.7 =0.8 =0.9 =1 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 u (t ) =0.5 =0.6 =0.7 =0.8 =0.9 =1 (b) Figure 16: (a) The numerical solutions of problem 4 using the FRK4 method with u0 = 0.5, h = 0.01 for different fractional orders β. (b) The numerical solutions of problem 4 using the FRK4 method with u0 = 0.5, h = 0.1 for different fractional orders β. Similar to previous cases, lower β values in this problem slow the convergence to equilibrium. Additionally, the smaller step size h in Figure (b) results in a smoother approximate solution. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 e r r o r a b s 10 -14 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 3 e r r o r a b s 10 -15 (b) Figure 17: (a) The errorabs of FRK4 for problem 4 with u0 = 0.5, h = 0.01 and orders β = 1. (b) The errorabs of FRK4 for problem 4 with u0 = 0.5, h = 0.0001 and orders β = 1. The errors are extremely low and are further reduced with h = 0.0001. The stable error profiles validate the precision of FRK4 for integer-order cases, establishing it as a reliable method for solving quadratic logistic models where high accuracy is essential. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 26 of 32 4.3. Problem 5 we consider the SEIR model defined in (2.4) with the parameter values given as b = 0.001555, p = 0.8, q = 0.95, r = 0.05, η = 0.05, θ = 0.002, γ = 0.003, d(N) = 0.00001+0.000007N (4.5) with the initial conditions S(0) = 140, E(0) = 0.01, I(0) = 0.02, N(0) = 141. (4.6) The given model is solved with h = 0.01, and final time t = 2500 we observe that the results presented in the Figures (18a-21) are in good agreement with the results in [42] 0 500 1000 1500 2000 2500 t 0 20 40 60 80 100 120 140 160 S u s p e n d e d =1 =0.95 =0.90 (a) 0 500 1000 1500 2000 2500 t 0 5 10 15 20 25 30 E x p o s e d =1 =0.95 =0.90 (b) Figure 18: (a) The plot shows the numerical solution for the susceptible population S for the problem defined in (2.4) with parameters values (4.5) and initial conditions given in (4.6). (b) The plot shows the numerical solution for the exposed population E for the problem defined in (2.4) with parameters values (4.5) and initial conditions given in (4.6). In Figure 18a, the susceptible class is expressed which has shown a sudden decrease in the population. The dynamics of the population is illustrated for three different fractional orders 1, 0.95 and 0.90. All these results have shown closer behavior in simulations with same pattern. A slight variation in the dynamics can easily be observed showing the impacts of the variation of the fractional order. For the least fractional order 0.9, we have observed that a more abrupt fall in the population has occur in the susceptible class in the early 500 units of time. This fall in the population is converted to the exposed population and can be seen in the Figure 18b. For the 0.90, the most early growth in the population of the class has been observed but has the same pattern of the dynamics as that for the order 1. While, after 250 to 500 units of time, the rise in the population is fallen down and have shifted to the infected population as shown in the Figure 19. There are variations in the dynamics of the infected class for the different fractional orders but have kept the Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 27 of 32 same pattern for all the three fractional orders. The stability in the population is observed after 500 for the order 0.90, after 750 for the order 0.95 and after 1000 units of time for the order 1. These three classes have finally grown to the recovered class where the same pattern in the population for the different orders can be observed in Figure 20. The total population of almost remains stable with a slight fall after infection which have increases the death tol in the population given in the Figure 21. 0 500 1000 1500 2000 2500 t 0 20 40 60 80 100 120 I n fe c ti o n s =1 =0.95 =0.90 Figure 19: The plot shows the numerical solution of infected population I for the problem defined in (2.4) with parameters values (4.5) and initial conditions given in (4.6). 0 500 1000 1500 2000 2500 t 0 20 40 60 80 100 120 R e c o v e r e d =1 =0.95 =0.90 Figure 20: The plot shows the numerical solution of the recovered population R for the problem defined in (2.4) with parameters values (4.5) and initial conditions given in (4.6). Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 28 of 32 0 500 1000 1500 2000 2500 t 10 30 50 70 90 110 130 150 170 N =1 =0.95 =0.90 Figure 21: The numerical solution N of the problem defined in (2.4) with parameters values (4.5) and initial conditions given in (4.6). 5. Conclusions In this article, we successfully developed and applied the FRK4 method to fractional differential equations (FDEs) arising in fractional logistic growth and SEIR models. The results demonstrate that the FRK4 method is a reliable and effective tool for approximat- ing solutions to FDEs in complex biological systems. By employing the Caputo derivative to describe fractional dynamics, we effectively captured the memory and hereditary char- acteristics inherent to such systems. The method proved particularly advantageous due to its high accuracy, flexibility in handling different types of fractional derivatives, and strong numerical stability, even in nonlinear contexts. However, these benefits come with an increased computational cost, especially when higher accuracy is sought through finer step sizes or more evaluations. This is primarily due to the non-local nature of fractional derivatives, which require integration over the entire history of the solution. Despite this computational overhead, our numerical experiments confirmed the ro- bustness and superior convergence properties of the FRK4 method compared to existing numerical approaches. This study not only validates the applicability of the FRK4 method to fractional logistic and SEIR models but also provides a solid foundation for its extension to a broader range of problems involving fractional dynamics. Acknowledgements The authors (H. Khan and J. Alzabut) would like to thank Prince Sultan University for their financial support to this work and its publication. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 29 of 32 Author Contributions K: Writing original draft, Supervision. S.M.: Validated the results, Conceptualization, Software. H.K.: Writing original draft, Investigation, Conceptualization. R.T.: Validated the results, Conceptualization, Software. J.A.: Supervision, Review manuscript, Editing. Declaration of Competing Interest The authors declare no competing interests. Data Availability The data is contained within the article. References [1] A Kilbas, Aleksandrovich, H M Srivastava, and J J Trujillo. Theory and applications of fractional differential equations, volume 204. elsevier, 2006. [2] I Podlubny. Fractional differential equations: an introduction to fractional deriva- tives, fractional differential equations, to methods of their solution and some of their applications, volume 198. elsevier, 1998. [3] C Drapaca. Fractional calculus in neuronal electromechanics. Journal of Mechanics of Materials and Structures, 12(1):35–55, 2016. [4] A Carpinteri and F Mainardi. Fractals and fractional calculus in continuum mechan- ics, volume 378. Springer, 2014. [5] L K B Kuroda, A V Gomes, R Tavoni, P F de Arruda Mancera, N Varalta, and R de Figueiredo Camargo. Unexpected behavior of caputo fractional derivative. Com- putational and Applied Mathematics, 36(3):1173–1183, 2017. [6] N Djeddi, S Hasan, M Al-Smadi, and S Momani. Modified analytical approach for generalized quadratic and cubic logistic models with caputo-fabrizio fractional deriva- tive. Alexandria Engineering Journal, 59(6):5111–5122, 2020. [7] S Arshad, A Sohail, and K Maqbool. Nonlinear shallow water waves: a fractional order approach. Alexandria Engineering Journal, 55(1):525–532, 2016. [8] X Liu, Kamran, and Y Yao. Numerical approximation of riccati fractional differential equation in the sense of caputo-type fractional derivative. Journal of Mathematics, 2020(1):1274251, 2020. [9] X J Yang, JAT Machado, and HM Srivastava. A new numerical technique for solving the local fractional diffusion equation: two-dimensional extended differential trans- form approach. Applied Mathematics and Computation, 274:143–151, 2016. Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 30 of 32 [10] X Yang, HM Srivastava, and C Cattani. Local fractional homotopy perturba- tion method for solving fractal partial differential equations arising in mathematical physics. Romanian Reports in Physics, 67(3):752–761, 2015. [11] F Mohammadi and C Cattani. A generalized fractional-order legendre wavelet tau method for solving fractional differential equations. Journal of Computational and Applied Mathematics, 339:306–316, 2018. [12] J Wang, Kamran, A Jamal, and X Li. Numerical solution of fractional-order fredholm integrodifferential equation in the sense of atangana–baleanu derivative. Mathematical Problems in Engineering, 2021(1):6662808, 2021. [13] X Qiang, Kamran, A Mahboob, and Y Chu. Numerical approximation of fractional-order volterra integrodifferential equation. Journal of Function Spaces, 2020(1):8875792, 2020. [14] S Momani and Z Odibat. Numerical comparison of methods for solving linear dif- ferential equations of fractional order. Chaos, Solitons & Fractals, 31(5):1248–1255, 2007. [15] S Arshad, A M Siddiqui, A Sohail, K Maqbool, and Z Li. Comparison of optimal homotopy analysis method and fractional homotopy analysis transform method for the dynamical analysis of fractional order optical solitons. Advances in Mechanical Engineering, 9(3):1687814017692946, 2017. [16] MHT Alshbool, AS Bataineh, I Hashim, and Osman Rasit Isik. Solution of fractional- order differential equations based on the operational matrices of new fractional bern- stein functions. Journal of King Saud University-Science, 29(1):1–18, 2017. [17] K Diethelm and G Walz. Numerical solution of fractional order differential equations by extrapolation. Numerical algorithms, 16:231–253, 1997. [18] K Diethelm and A D Freed. On the solution of nonlinear fractional-order differential equations used in the modeling of viscoplasticity. In Scientific computing in chemical engineering II: computational fluid dynamics, reaction engineering, and molecular properties, pages 217–224. Springer, 1999. [19] S Bhalekar and V Daftardar-Gejji. A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order. J. Fract. Calc. Appl, 1(5):1–9, 2011. [20] V Daftardar-Gejji, Y Sukale, and S Bhalekar. A new predictor–corrector method for fractional differential equations. Applied Mathematics and Computation, 244:158–182, 2014. [21] R Garrappa. On linear stability of predictor–corrector algorithms for fractional differ- ential equations. International Journal of Computer Mathematics, 87(10):2281–2290, 2010. [22] Z M Odibat and S Momani. An algorithm for the numerical solution of differen- tial equations of fractional order. Journal of Applied Mathematics & Informatics, 26(1 2):15–27, 2008. [23] P Rahimkhani, Y Ordokhani, and E Babolian. Numerical solution of fractional pan- tograph differential equations by using generalized fractional-order bernoulli wavelet. Journal of Computational and Applied Mathematics, 309:493–510, 2017. [24] Z Zhang, A Zeb, O Z Egbelowo, and V S Erturk. Dynamics of a fractional order Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 31 of 32 mathematical model for covid-19 epidemic. Advances in difference equations, 2020:1– 16, 2020. [25] M S Arshad, D Baleanu, M B Riaz, and M Abbas. A novel 2-stage fractional runge– kutta method for a time-fractional logistic growth model. Discrete Dynamics in Nature and Society, 2020(1):1020472, 2020. [26] SS Ezz-Eldien. On solving fractional logistic population models with applications. Computational and Applied Mathematics, 37:6392–6409, 2018. [27] S Das, PK Gupta, and K Vishal. Approximate approach to the das model of fractional logistic population growth. Applications and Applied Mathematics: An International Journal (AAM), 5(2):26, 2010. [28] S Z Rida, A A Farghaly, S A Azoz, and F Hussien. Global stability of a delayed fractional-order sei epidemic model with logistic growth. Appl. Math, 15(1):1–12, 2021. [29] R Abreu-Blaya, A Fleitas, J E Napoles Valdes, R Reyes, J M Rodriguez, and J M Sigarreta. On the conformable fractional logistic models. Mathematical Methods in the Applied Sciences, 43(7):4156–4167, 2020. [30] S M Ismail, L A Said, A A Rezk, A G Radwan, A H Madian, M F Abu-Elyazeed, and A M Soliman. Generalized fractional logistic map encryption system based on fpga. AEU-International Journal of Electronics and Communications, 80:114–126, 2017. [31] V V Tarasova and V E Tarasov. Logistic map with memory from economic model. Chaos, Solitons & Fractals, 95:84–91, 2017. [32] V Lakshmikantham and A S Vatsala. Theory of fractional differential inequalities and applications. Communications in Applied Analysis, 11(3-4):395–402, 2007. [33] E Demirci and N Ozalp. A method for solving differential equations of fractional order. Journal of computational and applied mathematics, 236(11):2754–2762, 2012. [34] C Milici, J Tenreiro Machado, and G Drăgănescu. Application of the euler and runge– kutta generalized methods for fde and symbolic packages in the analysis of some fractional attractors. International Journal of Nonlinear Sciences and Numerical Simulation, 21(2):159–170, 2020. [35] Constantin Milici, Gheorghe Drăgănescu, and J Tenreiro Machado. Introduction to fractional differential equations, volume 25. Springer, 2018. [36] Z A Khan, A Khan, T Abdeljawad, and H Khan. Computational analysis of fractional order imperfect testing infection disease model. Fractals, 30(05):2240169, 2022. [37] K Ramesh, A Khan, T Abdeljawad, et al. Analysis of the stability of a predator-prey model including the memory effect, double allee effect and holling type-i functional response. PloS one, 20(1):e0305179, 2025. [38] J Muhammad, U Younas, DK Almutairi, A Khan, and T Abdeljawad. Optical wave features and sensitivity analysis of a coupled fractional integrable system. Results in Physics, 68:108060, 2025. [39] R M Yaseen, N F Ali, Ahmed A Mohsen, A Khan, and T Abdeljawad. The modeling and mathematical analysis of the fractional-order of cholera disease: Dynamical and simulation. Partial Differential Equations in Applied Mathematics, 12:100978, 2024. [40] S Malik, K Muhammad, O Ahsan, M T Khan, R Sah, and Y Waheed. Advances in Kamran et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5944 32 of 32 zika virus vaccines and therapeutics: A systematic review. Asian Pacific Journal of Tropical Medicine, 17(3):97–109, 2024. [41] R Thinakaran, J Somasekar, V Neerugatti, and K Ganga. Advancements in skin can- cer detection: A comprehensive review of convolutional neural network approaches. In 2024 14th International Conference on Software Technology and Engineering (IC- STE), pages 232–236. IEEE, 2024. [42] Elif Demirci, Arzu Unal, et al. A fractional order seir model with density dependent death rate. Hacettepe journal of mathematics and statistics, 40(2):287–295, 2011. Introduction Preliminaries Fractional Runge-Kutta Method of Order Four Differential Equation of Fractional Order FRK4 Numerical Scheme Limitations and Challenges Fractional SEIR Model Results and Discussion Problem 1 Problem 2 Problem 3 Problem 4 Problem 5 Conclusions