EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 5684 ISSN 1307-5543 – ejpam.com Published by New York Business Global Computational Study of a Meshless Approach for Multi-Term Time-Fractional Models in Drug Dispersion and Absorption in Biological Tissues Imtiaz Ahmad1,∗, Rashid Jan2,3, Normy Norfiza Abdul Razak2, Aziz Khan4, Thabet Abdeljawad4,∗ 1 Institute of Informatics and Computing in Energy (IICE), Universiti Tenaga Nasional (UNITEN), Kajang 43000, Selangor, Malaysia 2 Institute of Energy Infrastructure (IEI), Department of Civil Engineering, College of Engineering, Universiti Tenaga Nasional (UNITEN), Putrajaya Campus, Jalan IKRAM-UNITEN, 43000 Kajang, Selangor, Malaysia 3 Mathematics Research Center, Near East University TRNC, Mersin 10, Nicosia 99138, Turkey 4 Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833, 11586 Riyadh, Saudi Arabia Abstract. This study presents a hybrid numerical approach for solving the multi-term time- fractional mobile-immobile diffusion equation. The proposed method combines a fractional Liouville- Caputo scheme for time derivatives with the embedded Runge–Kutta method (RK23) and employs a meshless technique using Fibonacci and Lucas polynomials for spatial derivatives. The accuracy and effectiveness of the approach are evaluated through numerical experiments using theMax error and RMS error norms. This study highlights the advantages of meshless techniques, particularly their flexibility in higher-dimensional applications and ease of implementation. The numerical re- sults obtained are validated against exact solutions, confirming the efficiency and accuracy of the proposed method. 2020 Mathematics Subject Classifications: 35R11, 65M70, 65L06, 76S05, 41A58 Key Words and Phrases: Hybrid numerical method, Mobile-immobile solute transport model, Caputo derivative, Embedded Runge–Kutta method, Fibonacci and Lucas polynomials 1. Introduction Fractional calculus is an advanced mathematical framework that generalizes the tra- ditional concepts of derivatives and integrals to fractional orders. This extension allows ∗Corresponding author. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5684 Email addresses: imtiazkakakhil@gmail.com (I. Ahmad), tabdeljawad@psu.edu.sa (T. Abdeljawad) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 2 of 24 for modeling systems with memory and hereditary characteristics, making it useful for processes that exhibit anomalous behaviors, such as diffusion, viscoelasticity, and biolog- ical systems. Unlike classical calculus, which limits its scope to integer-order operations, fractional calculus allows for more flexible and accurate modeling of complex real-world dynamics [1, 2]. As a result, it has garnered widespread applications across physics, engi- neering, finance, and signal processing [3–7]. Fractional differential equations (FDEs) expand the concept of differentiation to in- clude non-integer orders, thereby increasing the flexibility of classical differential equations [8]. This extension allows FDEs to capture a wide range of complex, memory-dependent processes in a variety of fields, such as mathematical biology, computational finance, and physics [9–11]. Fractional partial differential equations (PDEs) offer a robust framework for modeling phenomena with non-local dependencies and anomalous diffusion, which tra- ditional PDEs may inadequately capture [12]. In areas like fluid mechanics, for example, fractional PDEs have been effectively applied to simulate the behavior of viscoelastic and non-Newtonian fluids, where the presence of memory effects and non-linear stress-strain relationships demand a more complex approach than that offered by integer-order mod- els [13]. This adaptability makes FDEs a valuable tool for more realistic and accurate representations of complex, real-world systems. Anomalous diffusion, unlike classical diffusion with its even particle distribution over time, captures complex particle behaviors that lead in non-uniform spreading and extended- range interactions, especially useful for environments where standard models fall short. This complexity is modeled through the fractional diffusion equation, which incorpo- rates fractional derivatives to generalize the diffusion process, thus enabling the accu- rate simulation of varied and nonlocal particle movements observed in diverse systems, including porous media [12]. This approach is particularly advantageous in representing sub-diffusion in environments like groundwater or soil, where particles exhibit alternat- ing mobile and immobile states—a behavior captured by the mobile-immobile fractional diffusion equation variant [14, 15]. This model more accurately reflects the complex par- ticle dynamics in these settings by addressing the influences of diffusion, advection, and reaction processes on solute transport, resulting in a more comprehensive and realistic description of particle movement in heterogeneous systems. The concept of anomalous solute transport has important real-world applications, par- ticularly in environmental management and remediation. By providing insights into the irregular behavior of contaminants in groundwater, it facilitates the development of effec- tive remediation strategies. This model helps identify the sources of contamination, pre- dicting pollutant movement through soil and groundwater, and designing targeted cleanup methods [16]. Understanding subsurface water flow is essential for effective aquifer man- agement, as modeling can help optimize pumping and recharge rates, minimizing the risk of overuse and depletion [17]. In surface water systems, the model’s ability to predict helps in detecting pollution sources and developing strategies to reduce contamination. Additionally, it predicts contaminant and solvent migration in soils, contributing to ef- ficient land and water remediation efforts [18]. Anomalous solute transport modeling is also valuable in order to estimate carbon dioxide migration in underground formations, I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 3 of 24 offering insights into effective carbon sequestration methods to mitigate climate change impacts [19]. The mobile-immobile model (MIM) has become an essential framework for characteriz- ing non-Fickian solute transport in heterogeneous porous media, with important improve- ments from various researchers. Since the influential research work by Coats and Smith [20], the MIM has been applied to address challenges in hydrology, particularly in both unsaturated and saturated environments. Gao et al., [21] improved the MIM by incorpo- rating scale-dependent dispersion, effectively representing breakthrough curves in highly heterogeneous soil columns. They demonstrated that the liquid phase in porous media can be distinguish between mobile (flowing) and immobile (stagnant) regions based on their mobility, enabling a better representation of reactive solute transport. Additionally, Goltz and Roberts [22] derived three-dimensional solutions for solute transport in envi- ronments featuring both mobile and immobile zones, making their findings relevant for practical applications. Hasan et al., [23] used pore network modeling to investigate trans- port properties in two-phase flow, revealing that relative permeabilities could accurately determine stationary saturation. Zhang et al. [24] advanced this research by using lattice Boltzmann simulations to examine pore-scale solute transport across different types of porous media. Their findings indicated that materials with greater heterogeneity showed higher dispersion coefficients and exhibited non-Fickian transport behavior. They also highlighted the significant influence of mass transfer rates between mobile and stagnant regions on mechanical dispersion. The fractional mobile-immobile partial differential equation (PDE) models are increasingly applied in understanding drug dispersion and absorption in biological tissues, provide in- sights into the intricate phenomena associated with pharmaceutical delivery. These mod- els capture the non-Fickian behavior observed in biological systems, where drug molecules may exhibit varying degrees of mobility because of interactions with the tissue matrix and the existence of stagnant regions [25, 26]. For instance, they have been used to simu- late the transport of lipophilic drugs across skin layers, demonstrating that the fractional derivatives can accurately depict the delayed and anomalous diffusion processes that occur during permeation [27]. Additionally, these models aid in optimizing drug formulation and delivery systems by predicting the release profiles of drugs from carriers within biological tissues [28]. Furthermore, the fractional mobile-immobile approach has been used to study the kinetics of drug absorption in gastrointestinal tissues, yielding useful information for the design of oral drug delivery systems [29]. The adoption of fractional mobile-immobile PDE models enhances the understanding of drug transport mechanisms and makes it easier to create more effective therapeutic strategies. All these studies emphasize the importance of considering mobile and immobile regions in the analysis and modeling of solute transport within porous media. In this study, a two-term and a three-term time-fractional partial differential equa- tion are considered, and the proposed method is applied to the numerical solution of the I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 4 of 24 anomalous transport phenomenon in mobile-immobile systems [30, 31]. ∂P(X, T ) ∂T + µ ( ∂β1P(X, T ) ∂T β1 + ∂β2P(X, T ) ∂T β2 ) =α∆P(X, T ) + F (X, T ), T > 0, X = (X1, X2) ∈ Ω, 0 < β2 ≤ β1 ≤ 1, (1) and ∂P(X, T ) ∂T + µ ( ∂β1P(X, T ) ∂T β1 + ∂β2P(X, T ) ∂T β2 + ∂β3P(X, T ) ∂T β3 ) =α∆P(X, T ) + F (X, T ), T > 0, X = (X1, X2) ∈ Ω, 0 < β3 ≤ β2 ≤ β1 ≤ 1, (2) with conditions: P(X, 0) = 0, (3) P(X, T ) = g1(X, T ), X ∈ ∂Ω, (4) The Laplacian is denoted by ∆, with µ and α as established constants. Moreover, F (X, T ) and g1(X, T ) represent functional parameters. The Caputo fractional derivative [32] op- erator applied to the function P(X, T ) is as given below, and ∂βk ∂T βk , for k = 1, 2, 3. ∂βkP(X, T ) ∂T βk =  1 Γ(1−βk) T∫ 0 (T − η)−βkUη(X, η)dη, 0 < βk < 1 ∂P(X,T ) ∂T , βk = 1, in which Γ(.) is the gamma function. Partial differential equations (PDEs) are difficult to solve due to their natural com- plexity, which often limits the effectiveness of traditional solution methods. Consequently, there is an increasing need for efficient computational strategies that can deliver precise approximations. Among the many available techniques, meshless methods are distin- guished for their efficacy and accuracy, making them ideal for solving both fractional and non-fractional PDEs. Various meshless approaches exist, such as the reproducing kernel particle method, moving least squares method, natural element method, meshless radial basis function (RBF) methods, as well as methods based on Lucas and Fibonacci polynomials. The objective of this research is to investigate the underlying time-fractional models with a hybrid numerical approach. A meshless technique utilizing Lucas and Fibonacci polynomials is used to approximate spatial derivatives, while the Caputo definition com- bined with the embedded Runge–Kutta method is utilized to address the time fractional component. This method allows for the use of higher-order derivatives, improves accuracy even with fewer collocation points, and ultimately lowers computational costs. Further- more, it has been shown that these polynomials are very useful for resolving differential equations (DEs). Previous research emphasizes that effectively addressing boundary value I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 5 of 24 problems requires an understanding of the connections between Lucas and Chebyshev polynomials. For example, Lucas polynomials have been shown to be effective in solving higher-order DEs [33]. Furthermore, studies show that Volterra-Fredholm integral DEs can be effectively solved using Fibonacci polynomials [34]. Also, a hybrid approach that combines Taylor and Lucas polynomials has been proposed for solving delay difference equations [35]. The integration of hybrid Fibonacci and Lucas polynomial algorithms has led to novel solutions for time-dependent PDEs [36]. Combining Lucas polynomials with finite difference methods has produced efficient numerical schemes applicable to various PDE models [37–39]. 2. Basic Concepts in Polynomial Theory and Fractional Calculus Fractional calculus, which extends the concept of ordinary differentiation and inte- gration to non-integer orders, provides strong tools for simulating intricate, real-world phenomena with memory and hereditary properties [40–42]. Fractional derivatives have proven useful across diverse disciplines, including physics, biology, magnetohydrodynam- ics, engineering, and finance, due to their ability to capture dynamics. This framework is particularly useful in biological systems, such as drug dispersion and absorption in tissues, where the processes are often influenced by factors like irregular diffusion rates and varying absorption dynamics. Using fractional derivatives allows for a more precise representation of how drugs move through and interact with biological tissues, reflecting the fractal-like and memory-dependent nature of cellular environments. Definition 1: The Caputo’s fractional derivative [32]: ∂βkP(X, T ) ∂T βk = 1 Γ (1− βk) T∫ 0 ∂P(X, η) ∂η (T − η)−βk dη, 0 < βk < 1. (5) Definition 2: The Riemann-Liouville derivative [43]: ∂βkP(X, T ) ∂T βk = 1 Γ (1− βk) d dT T∫ T (P(X, ϱ)− P(X,T )) (ϱ− T )βk dϱ, 0 < βk < 1. (6) 2.1. Lucas and Fibonacci Polynomial Theory This section focuses on defining and utilizing Lucas and Fibonacci polynomials to ap- proximate unknown functions and their derivatives. Lucas polynomials [44]: Lucas polynomials can be expressed using a three-term recur- rence relation: Lj(X) = jLj−1(X) + Lj−2(X), j ≥ 2. (7) Assuming that L0(X) = 2 and L1(X) = X, Eq. (7) can generate a sequence of Lucas numbers with X = 1. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 6 of 24 Fibonacci polynomials [44]: Fibonacci polynomials, an extension of the Fibonacci num- bers, are defined by a recurrence relation involving three terms: Fj(X) = jFj−1(X) + Fj−2(X), j ≥ 2. (8) Beginning with F0(X) = 0 and F1(X) = 1, the sequence continues. For X = 1, Eq. (8) generates the classical Fibonacci number sequence. Lemma[44]: The mth derivative of the jth Lucas polynomial, denoted as Lj(X), can be expressed in terms of the jth Fibonacci polynomial, Fj(X), as follows: L (m) j (X) = jFj(X)Dm−1, Dm−1 = D ×D ×D · · · D︸ ︷︷ ︸ (m−1)time . (9) In this context, D denotes the (M + 1)× (M + 1) matrix, which is written as follows: D =  0 0 . . . 0 0 ... d 0  , where the computation of d follows [44]: djk = { j sin (k−j)π 2 , if k > j, 0, otherwise. 2.2. Approximation of Function Suppose that P(X) is a continuous function and that P ∈ L2(R). Under these condi- tions, P can be represented as a linear combination of the jth Lucas polynomials as given below: P(X) = ∞∑ j=0 ΛjLj(X). (10) Here, Λk represents the unknown coefficients, and Lj(X) denotes the Lucas polynomials. Similarly, P(X) can be expanded as a linear combination of the jth Fibonacci polynomials under the same conditions, as illustrated below: P(X) = ∞∑ j=0 ΛjFj(X). The unknown coefficients are represented by Λj , while the Fibonacci polynomials are written as Fj(X). I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 7 of 24 To compute the first-order derivative of the function P(X), it is expanded using the Lucas polynomial series. P ′(X) = ∞∑ j=0 ΛjL ′ j(X). (11) For the function P(X), the associated mth order derivative is as follows: Pm(X) = ∞∑ j=0 ΛjL (m) j (X), (12) where Pm(X) = dmP(X) dX m , L (m) j (X) = dmLj(X) dX m . The expressions in Eqs. (11) and (12) can be reformulated by applying the relationship provided in (9): P ′(X) = ∞∑ j=0 ΛjjFj(X), (13) Similarly, the formula below can be used to obtain the mth derivative: P(m)(X) = ∞∑ j=0 ΛjjFj(X)Dm−1, (14) in which D and Dm−1 are previously defined. Remark 1: Truncated series of Lucas and Fibonacci polynomials are frequently applied in numerical computations to represent P(X) along with its mth derivative. In particular, we consider the following: P(X) ≃ M∑ j=0 ΛjLj(X), M ∈ N, and P(m)(X) ≃ M∑ j=0 ΛjL (m) j (X) = M∑ j=0 ΛjjFj(X)Dm−1, M ∈ N. (15) 3. Proposed Method In this section, the proposed method for approximating the models presented in Eqs. (18)-(19) is discussed. For clarity and ease of notation within this section, we define the following: Pn+1(X) = P(X1, X2, T n+1), Pn+1 jk = P(X1j , X2k, T n+1), I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 8 of 24 where X1j , X2k denote the nodal points defined as below, and T n = n × δT , with δT signifying the size of the time step: X1j = a+ jhX1 , X2k = c+ khX2 , (j, k = 1, 2, · · · ,M, M ∈ N), In the space domain Ξ = [a, b]×[c, d] ⊂ R2, the mesh-size is represented by hX1 = (b−a)/M and hX2 = (d− c)/M in the X1 and X2 directions, respectively. 3.1. Temporal Discretization By using the established L1 formula, an approximation error of order O(δT 2−β1) is achieved in the discrete formulation at the (n+1)th time level, given that 0 < β1 ≤ 1 [38]. ∂β1P(X, T n+1) ∂T β1 = 1 Γ(1− β1) ∫ T n+1 0 ∂P(X, η) ∂η (T n+1 − η)−β1dη, = 1 Γ(1− β1) n∑ j=0 ∫ (j+1)×δT j×δT ∂P(X, η) ∂η (T n+1 − η)−β1dη, = 1 Γ(1− β1) n∑ j=0 [ Pj+1(X)− Pj(X) δT +O(δT ) ] ∫ (j+1)×δT j×δT ((j + 1)δT − η)−β1dη. Following integration, it gives: ∂β1P(X, T n+1) ∂T β1 =  Aβ1 ∑n j=0Kβ1(j) [ Pn−j+1(X)− Pn−j(X) ] +O(δT 2−β1), 0 < β1 < 1, Pn+1(X)−Pn(X) δT +O(δT ), β1 = 1, (16) where Aβ1 = δT −β1 Γ(2−β1) and Jβ1(j) = (j + 1)1−β1 − (j)1−β1 . Thus, by omitting the error term, we can write: ∂β1P(X, T n+1) ∂T β1 = Aβ1 [ Pn+1(X)− Pn(X) ] +Aβ1 n∑ j=1 Jβ1(j) [ Pn−j+1(X)− Pn−j(X) ] , (17) with Jβ1(j) = 1, j = 0. Before applying the θ-weighted scheme, let us express Eq. (1) by setting µ = 1: ∂P(X, T ) ∂T + ∂β1P(X, T ) ∂T β1 + ∂β2P(X, T ) ∂T β2 = LP(X, T ), X ∈ Ω, 0 < β2 ≤ β1 ≤ 1, T > 0, (18) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 9 of 24 with the conditions P(X, 0) = P0(X), P(X, T ) = f1(X, T ), X ∈ ∂Ω, (19) where L = α ( ∂2P(X1,X2,T ) ∂X2 1 + ∂P(X1,X2,T ) ∂X2 2 ) + F (X1, X2, T ). Next, we implement the θ-weighted scheme to transform Eq. (18) into: ∂P(X, T ) ∂T + ∂β1P(X, T ) ∂T β1 + ∂β2P(X, T ) ∂T β2 = θLPn+1(X) + (1− θ)LPn(X). (20) Using Eq. (16), we get: Pn+1(X)− Pn(X) δT +Aβ1 [ Pn+1(X)− Pn(X) ] +Aβ1 n∑ j=1 Jβ1(j) [ Pn−j+1(X)− Pn−j(X) ] +Aβ2 [ Pn+1(X)− Pn(X) ] +Aβ2 n∑ j=1 Jβ2(j) [ Pn−j+1(X)− Pn−j(X) ] = θLPn+1(X) + (1− θ)LPn(X). (21) Further simplification lead to P(n+1) =(I + δT Aβ1I + δT Aβ2I − δT θL)−1( (I + δT Aβ1I + δT Aβ2I + δT (1− θ)L)P(n)(X)−Qn β(X) ) , (22) where Qn β(X) = Aβ1 n∑ j=1 Jβ1(j) [ Pn−j+1(X)− Pn−j(X) ] +Aβ2 n∑ j=1 Jβ2(j) [ Pn−j+1(X)− Pn−j(X) ] . Thus far, Eq. (22) serves as the time-discrete representation of Eq. (18). Additionally, to achieve a accurate and computationally efficient solution for Eq. (22), the time direction is integrated with the Embedded Runge-Kutta method (RK23), which is outlined as follows [45, 46]: The Embedded Runge-Kutta method of orders 2 and 3 consists of Runge-Kutta meth- ods of these respective orders. To compute k1 and k2, we apply the second-order Runge- Kutta method (RK2), which provides a less precise solution at the (n + 1) time level, as shown below. P̂n+1 = Pn + δT 2 (k1 + k2). (23) The numerical solution P̂n+1 is compared to the RK3 solution Pn+1. If the difference between the two solutions exceeds a predefined threshold, the RK3 solution is discarded, and the process is restarted with a reduced step size ∆t. If the difference is sufficiently small, then V (n+1) is considered valid and accepted. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 10 of 24 The selection of the step size δT is based on the user-defined tolerance ϵ, as described below. e = ∥Pn+1 − P̂n+1∥ = ∥δT 3 (k1 − 2k3 + k2)∥. The local truncation errors associated with RK2 and RK3 methods can be expressed as P(1) exact − P̂(1) = C1δT 3 + . . . , P(1) exact − P(1) = C2δT 4 + . . . , where C1 and C2 are real constants. We can obtain from above equations: P(1) − P̂(1) = C1δT 3 + . . . , therefore e ≈ ∥C1∥δT 3. (24) To select a new step size, we use enew ≈ ∥C1∥δT 3 new < ϵ. (25) By comparing Eq. (25) with Eq. (24), we obtain ∥C1∥δT 3 new ∥C1∥δT 3 < ϵ e , which implies δTnew < δT ( ϵ e ) 1 3 . To fulfill the above inequality, we apply δT = 0.9δT ( ϵ e ) 1 3 . (26) The embedded RK23 method has a third-order accuracy of O(δT 3). 3.2. Space Discretization The 2D function P(X) = P(X1, X2) can be expressed as a truncated series of Lucas polynomials, represented as a product: P(X1, X2) ≃ M∑ j=0 M∑ m=0 ΛjmLj(X1)Lm(X2) = M∑ j=0 M∑ m=0 ΛjmLjm(X1, X2), (27) where Ljm(X1, X2) = Lj(X1)Lm(X2). Equation (27) can be written as: P(X1, X2) ≃ LT(X1, X2)Λ, (28) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 11 of 24 withΛ = [Λ00,Λ01,Λ10, · · · ,ΛMM ]T and LT(X1, X2) = [L00(X1, X2),L01(X1, X2),L10(X1, X2), · · · ,LMM (X1, X2)]. The partial derivatives of the function P(X1, X2) are determined as follows: ∂P(X1, X2) ∂X1 ≃ M∑ j=0 M∑ m=0 Λjm dLj(X1) dX1 Lm(X2) = M∑ j=0 M∑ m=0 Λjmj(Fj(X1))Lm(X2) = LT X1 (X1, X2)Λ, and ∂2P(X1, X2) ∂X2 1 ≃ M∑ j=0 M∑ m=0 Λjm d2Lj(X1) dX2 1 Lm(X2) = M∑ j=0 M∑ m=0 Λjmj(Fj(X1)D)Lm(X2) = LT X1X1 (X1, X2)Λ, where LT X1 (X1, X2) = { ∂ ∂X1 Ljm(X1, X2)}Mj,m=0 = {dLj(X1) dX1 Lm(X2)}Mj,m=0 = {j(Fj(X1))Lm(X2)}Mj,m=0, LT X1X1 (X1, X2) = { ∂2 ∂X2 1 Ljm(X1, X2)}Mj,m=0 = {d 2Lj(X1) dX2 1 Lm(X2)}Mj,m=0 = {j(Fj(X1)D)Lm(X2)}Mj,m=0. Similarly ∂P(X1, X2) ∂X2 ≃ M∑ j=0 M∑ m=0 ΛjmLj(X1) dLm(X2) dX2 = M∑ j=0 M∑ m=0 ΛjmmLjk(X1)(Fm(X2)) = LT X2 (X1, X2)Λ, and ∂2P(X1, X2) ∂X2 2 ≃ M∑ j=0 M∑ m=0 ΛjmLj(X1) d2Lm(X2) dX2 2 = M∑ j=0 M∑ m=0 ΛjmmLj(X1)(Fm(X2)D) = LT X2X2 (X1, X2)Λ, with LT X2 (X1, X2) = { ∂ ∂X2 Ljm(X1, X2)}Mj,m=0 = {Lj(X1) d2Lm(X2) dX2 2 }Mj,m=0 = {mLj(X1)(Fm(X2))}Mj,m=0. LT X2X2 (X1, X2) = { ∂2 ∂X2 2 Ljm(X1, X2)}Mj,m=0 = {Lj(X1) d2Lm(X2) dX2 2 }Mj,m=0 = {mLj(X1)(Fm(X2)D)}Mj,m=0. These results can be expanded upon in the following ways: ∆P(X1, X2) = ∂2P ∂X2 1 + ∂2P ∂X2 2 ≃ ( LT X1X1 (X1, X2) + LT X2X2 (X1, X2) ) Λ. (29) Now, time-dependent function P(X1, X2, T ) and its partial derivatives can be express as: P(X1, X2, T ) ≈ LT(X1, X2)Λ(T ), (30) I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 12 of 24 ∂2P(X1, X2, T ) ∂X2 1 ≈ LT X1X1 (X1, X2)Λ(T ), (31) ∂2P(X1, X2, T ) ∂X2 2 ≈ LT X2X2 (X1, X2)Λ(T ), (32) with ∆P(X1, X2, T ) ≈ ( LT X1X1 (X1, X2) + LT X2X2 (X1, X2) ) Λ(T ). (33) Here, the unknown coefficients are denoted by Λ(T ) of the vector of time-dependent. By combining the information obtained from Eq. (33) with Eq. (22), along with the appropriate boundary condition, the necessary numerical results can be derived. 4. Analysis of Convergence and Stability For the purpose of error analysis, we will use the theorem outlined below. Theorem: Let P̂ denote the exact solution and P the computed solution obtained using the pro- posed meshless method for the underlying problem. Then, the error |E| can be expressed as: |E| = |P̂− P|, and is bounded by [38]: |E| ≤ 4e(2λ) cosh2(2R)λ2(M+1) ((M + 1)!)2 . Proof: Let the absolute error of the proposed method is: |E| = |P̂− P|, where P̂ = ∑∞ j=0 ∑∞ m=0 ΛjΛmLj(X1)Lm(X2) and P = ∑M j=0 ∑M m=0 ΛjΛmLj(X1)Lm(X2). Next, the term to be truncated is: |E| = ∞∑ j=M+1 ∞∑ m=M+1 ΛjΛmLj(X1)Lm(X2) (34) As shown in [47], that: Lm(r) ≤ 2ϱm, |Λm| ≤ Rm cosh(2R) m! , I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 13 of 24 where the well-known golden ratio is represented by ϱ. Thus, from Eq. (34), it follows that: |E| = 4 cosh2(2R) ∞∑ j=M+1 ∞∑ m=M+1 λj+m j!m! , where λ = Pϱ. Which can be written as: |E| = 4e(2λ) cosh2(2R) [ 1− Γ(M + 1, λ) Γ(M + 1) ]2 . (35) In this context, the incomplete gamma function is denoted as Γ(M + 1, λ), while the complete gamma function is represented as Γ(M + 1) [48]. The integral form of Equation (35) is given by: |E| ≤ 4e(2λ) cosh2(2R) (M !)2 [∫ λ 0 T Me(−T )dT ]2 . Since e−T < 1, therefore, ∀T > 0, we obtain: |E| ≤ 4e(2λ) cosh2(2R)λ2(M+1) ((M + 1)!)2 . 4.1. Stability analysis Stability analysis in this study is performed using a matrix technique, following the approach outlined in earlier research [49, 50]. Theorem: Let P represent the approximated solution obtained for the given problem. In this scenario, 𭟋 = LG−1HL−1 is the amplification matrix. Stability is guaranteed if the maximum absolute eigenvalue (ρ(𭟋)) of 𭟋 is less than or equal to 1. Proof: To demonstrate that 𭟋 = LG−1HL−1, we start by expressing from Eq. (15): Λn = L−1Pn. (36) The relationship can be expressed as follows: Pn+1 = LG−1HL−1Pn + LG−1Qn+1. (37) Furthermore, the error vector: E = P̂− P, (38) will satisfy the following condition: En+1 = 𭟋En. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 14 of 24 The amplification matrix is given by 𭟋 = LG−1HL−1. For the method to remain stable, it must satisfy the Lax–Richtmyer stability criterion [51], which states that: ||𭟋|| ≤ 1. 5. Numerical Results This section demonstrates the accuracy and effectiveness of the proposed technique for computing the solution of the multi-term time-fractional mobile-immobile diffusion equation. The method uses Lucas and Fibonacci polynomials as basis functions. Unless stated otherwise, three test problems are used to assess the performance of the method. The time step is set as δT = 0.0005, and the spatial domain is [0, 1]. The accuracy of the method is determined by the following formula: Absolute− error = |P̂− P|, Max error = max(Absolute− error), RMS = √√√√∑Nn i=1 ( P̂i − Pi )2 N . (39) Problem 1. The models described by equations (1)-(2), where µ = α = 1, have an exact solution given by: P(X, T ) = T 2 cos(2πX1) cos(2πX2), X = (X1, X2) ∈ Ω, (40) The numerical results for Problem 1 are generated using the proposed method and are presented in Table 1. The final time T = 0.5, number of collocation points M , time step δT and fractional orders β1 = β2 = 0.5 for the two-term case and β1 = β2 = β3 = 0.5 for the three-term case are considered. Also, the accuracy is measured through the Max error and RMS error norms. The results presented in this table show the superior performance of the proposed method. The simulation results for various T values, δT and M = 20 are presented in Table 2. The outcomes of this scenario indicate a noticeable improvement in accuracy. A comparison of the proposed method with the meshless RBF technique [30] for different time-fractional orders is shown in Table 3. As shown in the table, the proposed method demonstrates superior accuracy compared to the one in [30]. The comparison between the exact and numerical solutions for the two-term model, along with the absolute error norm, is shown in Figure 1, highlighting the method’s commendable accuracy. Similarly, the contour plots for the two-term and three-term models are presented in Figure 2. Problem 2. The models described by equations (1)-(2), where µ = α = 1, have an exact solution given by: P(X, T ) = T 2e(X1+X2), X = (X1, X2) ∈ Ω, (41) The numerical results for T = 0.5 are presented for various fractional orders βk and M , as detailed in Table 4. For the case of two terms, the fractional order is set as β = β1 = β2, I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 15 of 24 Table 1: Problem 1, Results for different numbers of nodal points and time-step sizes. M = 10 M = 20 M = 30 δT RMS Max error RMS Max error RMS Max error Two-term 0.05 9.7392e-03 6.8256e-02 8.4581e-03 6.2623e-02 6.5287e-03 5.3262e-02 β1 = β2 = 0.5 0.005 1.8351e-03 6.5328e-03 3.5518e-03 7.1735e-03 3.0124e-03 6.6805e-03 0.0005 6.7391e-04 7.3955e-04 6.2784e-04 7.0381e-04 5.5531e-04 7.1839e-04 Three-term 0.05 9.1613e-03 6.0534e-02 6.9831e-03 5.7504e-02 5.3843e-03 4.1438e-02 β1 = β2 = β3 = 0.5 0.005 1.0484e-03 6.0514e-03 2.5388e-03 5.0379e-03 1.8459e-03 5.5633e-03 0.0005 5.2719e-04 6.3736e-04 5.2857e-04 6.8145e-04 4.0015e-04 4.9745e-04 Table 2: Problem 1, Results for different time-step sizes and final time. T = 0.5 T = 1 δT RMS Max error RMS Max error Two-term 0.025 6.2463e-03 4.5347e-02 6.1248e-02 1.3648e-01 β1 = 0.75, β2 = 0.5 0.0025 1.4304e-03 5.4592e-03 6.2048e-03 3.6782e-02 0.00025 5.7392e-04 7.0472e-04 8.2943e-04 6.1320e-03 Three-term 0.025 5.6935e-03 3.6718e-02 5.5789e-02 8.5726e-02 β1 = 0.75, β2 = 0.5, β3 = 0.25 0.0025 1.1075e-03 4.0431e-03 3.4934e-03 2.2903e-02 0.00025 4.1042e-04 5.9832e-04 6.47253e-04 4.0237e-03 −1 −0.5 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 X1 = X2 S o lu t io n Exact Appro −1 −0.5 0 0.5 1 0 0.5 1 1.5 2 2.5 3 x 10 −3 X1 = X2 A b s o lu t e e r r o r n Figure 1: Problem 1, a comparison between the the numerical and exact solutions for the two-term model equation, along with the absolute error norm. while for three terms, β = β1 = β2 = β3. The table clearly indicates that accuracy improves with an increase in the number of nodes and the terms in the time-fractional order.Additionally, Figures 3–4 illustrate comparisons between the two-term and three- I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 16 of 24 Table 3: Problem 1, a comparison of the maximum error of the proposed method with the meshless RBF method [30] is presented for various fractional orders. Two-term Three-term βk proposed method meshless RBF method [30] proposed method meshless RBF method [30] 0.25 6.2849e-04 3.1321e-03 5.6248e-04 3.1176e-03 0.5 6.1355e-04 3.0031e-03 4.1736e-04 2.9294e-03 0.75 2.1426e-04 2.3804e-03 1.4510e-04 2.1178e-03 0.9 9.9826e-05 1.6083e-03 8.0655e-05 1.2919e-03 0 .0 0 0 4 5 4 3 8 0 .0 0 0 4 5 4 3 8 0.00045438 0 .0 0 0 4 5 4 3 8 0.00045438 0.0 0045438 0.00045438 0.00045438 0 .0 0 0 4 5 4 3 8 0.00045438 0.00045438 0.00045438 0.00045438 0.00045438 0 .0 0 0 4 5 4 3 8 0.00045438 0 .0 0 0 4 5 4 3 8 0.00045438 0 .0 0 0 4 5 4 3 8 0.00045438 0.00045438 0 .0 0 0 4 5 4 3 8 0.00090875 0.00090875 0 .0 0 0 9 0 8 7 5 0.00090875 0.0 0090875 0 .0 0 0 9 0 8 7 5 0 .0 0 0 9 0 8 7 5 0.00090875 0.00090875 0 .0 0 0 9 0 8 7 5 0 .0 0 0 9 0 8 7 5 0.00090875 0 .0 0 0 9 0 8 7 5 0.00090875 0.0013631 0 .0 0 1 3 6 3 1 0 .0 0 1 3 6 3 1 0.0013631 0 .0 0 1 3 6 3 1 0.0013631 0.0013631 0.0 013631 0.0013631 0.0013631 0.0013631 0.0013631 0.0018175 0 .0 0 1 8 1 7 5 0.0018175 0 .0 0 1 8 1 7 5 0 .0 0 1 8 1 7 5 0 .0 0 1 8 1 7 5 0.0018175 0.0018175 0.0018175 0.0018175 0 .0 0 2 2 7 1 9 0 .0 0 2 2 7 1 9 0.0022719 0.0022719 X1 X 2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x 10 −3 0 .0 0 0 4 4 2 9 8 0 .0 0 0 4 4 2 9 8 0 .0 0 0 4 4 2 9 8 0.0 0044298 0 .0 0 0 4 4 2 9 8 0.00044298 0.00044298 0.00044298 0.00044298 0.00044298 0 .0 0 0 4 4 2 9 8 0.00044298 0.00044298 0 .0 0 0 4 4 2 9 8 0.00044298 0 .0 0 0 4 4 2 9 8 0.00044298 0 .0 0 0 4 4 2 9 8 0.00044298 0.00044298 0.00044298 0.00044298 0.00088596 0.00088596 0.0 0088596 0.00088596 0.0 0088596 0 .0 0 0 8 8 5 9 6 0 .0 0 0 8 8 5 9 6 0.00088596 0.00088596 0 .0 0 0 8 8 5 9 60 .0 0 0 8 8 5 9 6 0.00088596 0 .0 0 0 8 8 5 9 6 0.00088596 0.0013289 0 .0 0 1 3 2 8 9 0 .0 0 1 3 2 8 9 0.0013289 0 .0 0 1 3 2 8 9 0.0013289 0.0013289 0.0 013289 0.0013289 0.0013289 0.0013289 0.0013289 0.0017719 0 .0 0 1 7 7 1 9 0.0017719 0 .0 0 1 7 7 1 9 0 .0 0 1 7 7 1 9 0 .0 0 1 7 7 1 9 0.0017719 0.0017719 0 .0 0 2 2 1 4 9 0 .0 0 2 2 1 4 9 0.0022149 0.0022149 X1 X 2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x 10 −3 Figure 2: Problem 1, the absolute error norm for the two-term and three-term model equations, as obtained using the proposed method. term examples. These figures confirm the accuracy of the proposed approach by displaying both the numerical and exact solutions, along with the corresponding absolute errors. To evaluate the efficiency of the proposed method on non-rectangular domains with uniform and non-uniform data points, the results for Test Problem 2 are presented in Figures 5–6. These findings demonstrate that, regardless of the shape of the computational domain, the method achieves appropriate accuracy. Problem 3. The models described by equations (1)-(2), where µ = α = 1, have an exact solution given by: P(X, T ) = T 2 sin(2πX1) sin(2πX2), X = (X1, X2) ∈ Ω, (42) Table 5 shows the results of suggested method for Problem 3. A comparison is made between the proposed method and the one outlined in [30] for various fractional orders. The final time T = 0.5 is considered to be fixed. For the two-term case, β = β1 = β2 is used, and for the three-term case, β = β1 = β2 = β3. The results demonstrate that the proposed method provides more accurate outcomes compared to the RBF method presented in [30]. The table also indicates that more accurate numerical results are achieved by increasing the number of nodes and terms in the time-fractional derivative. Figure 7 presents a contour I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 17 of 24 Table 4: Problem 1, Results for different nodal points and fractional-orders. Two-term Three-term M RMS Max error RMS Max error β = 0.25 10 6.8537e-04 7.1453e-04 5.8534e-04 6.8673e-04 20 5.7357e-05 5.4275e-04 4.5236e-05 4.4433e-04 30 4.9364e-05 4.7544e-04 4.0249e-05 3.3375e-04 β = 0.5 10 6.6063e-04 7.4643e-04 5.6296e-04 6.4126e-04 20 5.0953e-05 5.6489e-04 4.3343e-05 4.6455e-04 30 3.5467e-05 4.2314e-04 2.5979e-05 3.3652e-04 β = 0.75 10 6.1464e-04 6.6234e-04 5.7264e-04 5.8596e-04 20 5.0116e-05 4.9712e-04 4.0234e-05 3.2297e-04 30 3.0277e-05 3.3556e-04 2.9634e-05 1.3538e-04 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.5 1 1.5 2 X1X2 S o lu t io n Exact Appro −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.5 1 1.5 2 2.5 3 x 10 −3 X1 X2 A b s o lu t e e r r o r Figure 3: Problem 2, a comparison between the the numerical and exact solutions for the two-term model equation, along with the absolute error norm. map comparing the two-term and three-term cases, while Figure 8 shows a comparison between the numerical and exact solutions. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 18 of 24 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.5 1 1.5 2 X1X2 S o lu t io n Exact Appro −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.5 1 1.5 2 2.5 x 10 −3 X1 X2 A b s o lu t e e r r o r Figure 4: Problem 2, a comparison between the numerical and exact solutions for the three-term model equation, along with the absolute error norm. 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X1 X 2 0.00027325 0 .0 0 0 2 7 3 2 5 0.00027325 0.00027325 0 .0 0 0 2 7 3 2 5 0.00027325 0 .0 0 0 2 7 3 2 5 0.00027325 0.0 0027325 0.00027325 0 .0 0 0 2 7 3 2 5 0.00027325 0.00027325 0 .0 0 0 2 7 3 2 5 0 .0 0 0 2 7 3 2 5 0.0005465 0 .0 0 0 5 4 6 5 0 .0 0 0 5 4 6 5 0 .0 0 0 5 4 6 5 0.0005465 0 .0 0 0 5 4 6 5 0.0005465 0 .0 0 0 5 4 6 5 0.0005465 0 .0 0 0 5 4 6 5 0 .0 0 0 5 4 6 5 0.0005465 0.0005465 0.00081975 0.0 0081975 0 .0 0 0 8 1 9 7 5 0 .0 0 0 8 1 9 7 5 0.00081975 0.00081975 0.00081975 0.00081975 0.00081975 0 .0 0 0 8 1 9 7 5 0.00081975 0.001093 0 .0 0 1 0 9 3 0 .0 0 1 0 9 3 0 .0 0 1 0 9 3 0 .0 0 1 0 9 3 0.001093 0.001093 0.001093 0.001093 0.0013662 0 .0 0 1 3 6 6 2 0 .0 0 1 3 6 6 20 .0 0 1 3 6 6 2 0.0013662 0.0013662 0.0013662 0.0013662 0 .0 0 1 6 3 9 5 0 .0 0 1 6 3 9 5 0.0016395 0.0016395 0.0016395 0.0016395 0 .0 0 1 9 1 2 7 0.0019127 0.0019127 0.0019127 0 .0 0 2 1 8 6 0 .0 0 2 1 8 6 0.002186 0.002186 X1 X 2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 −3 Figure 5: Problem 2, computational domain with absolute error norm. −1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 X1 X 2 0 .0 0 0 1 7 2 1 7 0 .0 0 0 1 7 2 1 7 0.00017217 0.00017217 0.00017217 0 .0 0 0 1 7 2 1 7 0.00017217 0.00017217 0.0 0034435 0.0 00344350.00034435 0 .0 0 0 3 4 4 3 5 0 .0 0 0 3 4 4 3 5 0.0 0034435 0 .0 0 0 3 4 4 3 5 0.0 0051652 0 .0 0 0 5 1 6 5 2 0.00051652 0.00051652 0 .0 0 0 5 1 6 5 2 0.00051652 0.00051652 0.00051652 0.0006887 0.0006887 0 .0 0 0 6 8 8 7 0.0006887 0.0006887 0 .0 0 0 6 8 8 7 0.0006887 0.0006887 0.00086087 0.00086087 0.0 0086087 0.00086087 0.001033 0.001033 0.001033 0.001033 0.0012052 0.0012052 0 .0 0 1 2 0 5 2 0.0 013774 X1 X 2 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 x 10 −4 Figure 6: Problem 2, computational domain with absolute error norm. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 19 of 24 Table 5: Problem 3, a comparison of the maximum error of the proposed method with the meshless RBF method [30] is presented for different fractional orders and nodal points. Two-term Three-term M proposed method meshless RBF method [30] proposed method meshless RBF method [30] β = 0.25 10 5.4431e-04 2.8221e-03 4.7425e-04 2.7726e-03 20 7.1847e-05 7.6845e-04 6.4253e-05 7.5473e-04 30 6.3638e-05 3.5439e-04 5.1349e-05 3.4798e-04 β = 0.45 10 5.5394e-04 2.7834e-03 5.1436e-04 2.7166e-03 20 6.5637e-05 7.5756e-04 5.5207e-05 7.3900e-04 30 5.1328e-05 3.4906e-04 4.7508e-05 3.4028e-04 β = 0.65 10 5.2443e-04 2.7262e-03 4.1438e-04 2.6357e-03 20 6.0675e-05 7.4106e-04 5.5501e-05 7.1584e-04 30 4.9752e-05 3.4031e-04 3.1683e-05 3.2813e-04 X1 X 2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 x 10 −4 X1 X 2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 x 10 −4 Figure 7: Problem 3, the absolute error norm for the two-term and three-term model equations, as obtained using the proposed method. −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −0.4 −0.2 0 0.2 0.4 X1X2 S o lu t io n Exact Appro Figure 8: Problem 3, a comparison between the numerical and exact solutions for the two-term model equation. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 20 of 24 6. Conclusions In this study, an efficient numerical scheme is utilized for the multi-term time-fractional mobile-immobile diffusion equation, which effectively characterizes anomalous solute trans- port processes. Our approach integrates the fractional Liouville-Caputo scheme for han- dling time derivatives with the embedded Runge–Kutta method (RK23), leveraging a meshless framework based on Lucas and Fibonacci polynomials for spatial derivatives. The semi-discretized problem is addressed through a hybrid numerical approach using Lu- cas and Fibonacci polynomials, enabling high-accuracy approximations of the governing equations. The method is tested through successful applications to three test problems, with comprehensive tables and figures showcasing the numerical results, further demon- strating its effectiveness and emphasizing its potential for real-world applications in simu- lating anomalous transport phenomena within complex systems. This research highlights the advantages of the hybrid approach, particularly its flexibility in higher-dimensional scenarios, its straightforward implementation, and its applicability to real-world problems in diverse domains such as hydrology, bioengineering, and environmental modeling. Acknowledgements and Funding Aziz Khan and Thabet Abdeljawad would like to thank Prince Sultan University for paying the APC and supporting through TAS research lab. Declarations Competing interests: Authors declare no conflict of interest. Availability of data and materials: Data will be provided on request to the corresponding author. References [1] Fulgence Mansal and Ndolane Sene. Analysis of fractional fishery model with reserve area in the context of time-fractional order derivative. Chaos, Solitons & Fractals, 140:110200, 2020. [2] Ndolane Sene and José Francisco Gómez Aguilar. Fractional mass-spring-damper system described by generalized fractional order derivatives. fractal and fractional, 3(3):39, 2019. [3] Imtiaz Ahmad, Asmidar Abu Bakar, Hijaz Ahmad, Aziz Khan, and Thabet Abdel- jawad. Investigating virus spread analysis in computer networks with Atangana– Baleanu fractional derivative models. Fractals, 2440043:17, 2024. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 21 of 24 [4] Imtiaz Ahmad, Ihteram Ali, Rashid Jan, Sahar Ahmed Idris, and Mohamed Mousa. Solutions of a three-dimensional multi-term fractional anomalous solute transport model for contamination in groundwater. Plos one, 18(12):e0294348, 2023. [5] Ziad Ur Rehman, Salah Boulaaras, Rashid Jan, Imtiaz Ahmad, and Salma Bahra- mand. Computational analysis of financial system through non-integer derivative. Journal of Computational Science, 75:102204, 2024. [6] Imtiaz Ahmad, Ibrahim Mekawy, Muhammad Nawaz Khan, Rashid Jan, and Salah Boulaaras. Modeling anomalous transport in fractal porous media: A study of fractional diffusion PDEs using numerical method. Nonlinear Engineering, 13(1):20220366, 2024. [7] Salah Boulaaras, Rashid Jan, Amin Khan, Ali Allahem, Imtiaz Ahmad, and Salma Bahramand. Modeling the dynamical behavior of the interaction of T-cells and human immunodeficiency virus with saturated incidence. Communications in Theoretical Physics, 76(3):035001, 2024. [8] Anatolĭı Kilbas. Theory and applications of fractional differential equations. [9] RL Magin. Fractional calculus in bioengineering begell house publishers. Inc., Con- necticut, 2006. [10] Mustafa Inc, Muhammad Nawaz Khan, Imtiaz Ahmad, Shao-Wen Yao, Hijaz Ahmad, and Phatiphat Thounthong. Analysing time-fractional exotic options via efficient local meshless method. Results in Physics, 19:103385, 2020. [11] Hijaz Ahmad, Muhammad Nawaz Khan, Imtiaz Ahmad, Mohamed Omri, and Maged F Alotaibi. A meshless method for numerical solutions of linear and non- linear time-fractional Black-Scholes models. AIMS Math, 8(8):19677–19698, 2023. [12] Ralf Metzler and Joseph Klafter. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics reports, 339(1):1–77, 2000. [13] Michele Caputo and Mauro Fabrizio. A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation & Applications, 1(2):73–85, 2015. [14] Brian Berkowitz, Andrea Cortis, Marco Dentz, and Harvey Scher. Modeling non- Fickian transport in geological formations as a continuous time random walk. Reviews of Geophysics, 44(2), 2006. [15] Imtiaz Ahmad, Muhammad N Khan, Mustafa Inc, Hijaz Ahmad, and KS Nisar. Nu- merical simulation of simulate an anomalous solute transport model via local meshless method. Alexandria Engineering Journal, 59:2827–2838, 2020. [16] Brian Berkowitz, Harvey Scher, and Stephen E Silliman. Anomalous transport in laboratory-scale, heterogeneous porous media. Water Resources Research, 36(1):149– 158, 2000. [17] Lynn W Gelhar, Claire Welty, and Kenneth R Rehfeldt. A critical review of data on field-scale dispersion in aquifers. Water resources research, 28(7):1955–1974, 1992. [18] Chunmiao Zheng, Gordon D Bennett, et al. Applied contaminant transport modeling, volume 2. Wiley-Interscience New York, 2002. [19] Ruben Juanes, EJ Spiteri, FM Orr Jr, and MJ Blunt. Impact of relative permeability hysteresis on geological CO2 storage. Water resources research, 42(12), 2006. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 22 of 24 [20] KH Coats and BD Smith. Dead-end pore volume and dispersion in porous media. Society of petroleum engineers journal, 4(01):73–84, 1964. [21] Guangyao Gao, Hongbin Zhan, Shaoyuan Feng, Bojie Fu, Ying Ma, and Guan- hua Huang. A new mobile-immobile model for reactive solute transport with scale- dependent dispersion. Water Resources Research, 46(8), 2010. [22] Mark N Goltz and Paul V Roberts. Three-dimensional solutions for solute transport in an infinite medium with mobile and immobile zones. Water Resources Research, 22(7):1139–1148, 1986. [23] Sharul Hasan, Vahid Joekar-Niasar, Nikolaos K Karadimitriou, and Muhammad Sahimi. Saturation dependence of non-Fickian transport in porous media. Water Resources Research, 55(2):1153–1166, 2019. [24] Chunwei Zhang, Kazuki Kaito, Yingxue Hu, Anindityo Patmonoaji, Shintaro Mat- sushita, and Tetsuya Suekane. Influence of stagnant zones on solute transport in heterogeneous porous media at the pore scale. Physics of Fluids, 33(3), 2021. [25] Richard L Magin. Fractional calculus models of complex dynamics in biological tis- sues. Computers & Mathematics with Applications, 59(5):1586–1593, 2010. [26] Vesna Miskovic-Stankovic and Teodor M Atanackovic. On a system of equations with general fractional derivatives arising in diffusion theory. Fractal and Fractional, 7(7):518, 2023. [27] Michele Caputo and Cesare Cametti. Fractional derivatives in the transport of drugs across biological materials and human skin. Physica A: Statistical Mechanics and its Applications, 462:705–713, 2016. [28] Sanxiu Lu, W Fred Ramirez, and Kristi S Anseth. Modeling and optimization of drug release from laminated polymer matrix devices. AIChE journal, 44(7):1689– 1696, 1998. [29] Po-Chang Chiang, Jia Liu, Peter Fan, and Harvey Wong. Exploring a kinetic model approach in biopharmaceutics: Estimating the fraction absorbed of orally adminis- tered drugs in humans. Journal of Pharmaceutical Sciences, 107(7):1798–1805, 2018. [30] Hanaa Abu-Zinadah, MD Alsulami, and Hijaz Ahmad. Application of efficient hybrid local meshless method for the numerical simulation of time-fractional PDEs arising in mathematical physics and finance. The European Physical Journal Special Topics, 232(14):2595–2605, 2023. [31] Q. Liu, F. Liu, I. Turner, V. Anh, and Y. T. Gu. A RBF meshless approach for mod- eling a fractal mobile/immobile transport model. Applied Mathematics and Compu- tation, 226:336–347, 2014. [32] Michele Caputo. Linear models of dissipation whose Q is almost frequency independent-II. Geophysical Journal International, 13(5):529—-539, 1967. [33] Muhammed Çetin, Mehmet Sezer, and Coşkun Güler. Lucas polynomial approach for system of high-order linear differential equations and residual error estimation. Mathematical problems in engineering, 2015, 2015. [34] Farshid Mirzaee and Seyede Fatemeh Hoseini. Application of Fibonacci collocation method for solving Volterra–Fredholm integral equations. Applied Mathematics and Computation, 273:637–644, 2016. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 23 of 24 [35] N Baykuş-Savaşaneril and Mehmet Sezer. Hybrid Taylor-Lucas collocation method for numerical solution of high-order Pantograph type delay differential equations with variables delays. Applied Mathematics and Information Sciences, 11(6):1795–1801, 2017. [36] Ömer Oruç. A new numerical treatment based on Lucas polynomials for 1D and 2D sinh-Gordon equation. Communications in Nonlinear Science and Numerical Simulation, 57:14–25, 2018. [37] Imtiaz Ahmad, Sakhi Zaman, et al. Local meshless differential quadrature collocation method for time-fractional PDEs. Discrete & Continuous Dynamical Systems-Series S, 13(10), 2020. [38] Imtiaz Ahmad, Abdulrahman Obaid Alshammari, Rashid Jan, Normy Norfiza Ab- dul Razak, and Sahar Ahmed Idris. An efficient numerical solution of a multi- dimensional two-term fractional order PDE via a hybrid methodology: The Caputo– Lucas–Fibonacci approach with strang splitting. Fractal and Fractional, 8(6):364, 2024. [39] Imtiaz Ahmad, Asmidar Abu Bakar, Ihteram Ali, Sirajul Haq, Salman Yussof, and Ali Hasan Ali. Computational analysis of time-fractional models in energy infrastruc- ture applications. Alexandria Engineering Journal, 82:426–436, 2023. [40] Hasib Khan, Altaf Hussain Rajpar, Jehad Alzabut, Muhammad Aslam, Sina Etemad, and Shahram Rezapour. On a fractal–fractional-based modeling for influenza and its analytical results. Qualitative Theory of Dynamical Systems, 23(2):70, 2024. [41] Saim Ahmed, Ahmad Taher Azar, Mahmoud Abdel-Aty, Hasib Khan, and Jehad Alz- abut. A nonlinear system of hybrid fractional differential equations with application to fixed time sliding mode control for Leukemia therapy. Ain Shams Engineering Journal, 15(4):102566, 2024. [42] Hasib Khan, Jehad Alzabut, JF Gómez-Aguilar, and Abdulwasea Alkhazan. Essential criteria for existence of solution of a modified-ABC fractional order smoking model. Ain Shams Engineering Journal, 15(5):102646, 2024. [43] Guy Jumarie. Stock exchange fractional dynamics defined as fractional exponential growth driven by (usual) gaussian white noise. application to fractional Black-Scholes equations. Insurance: Mathematics and Economics, 42(1):271–287, 2008. [44] Ihteram Ali, Sirajul Haq, Kottakkaran Sooppy Nisar, and Dumitru Baleanu. An effi- cient numerical scheme based on Lucas polynomials for the study of multidimensional Burgers-type equations. Advances in Difference Equations, 2021(1):1–24, 2021. [45] Imtiaz Ahmad, Siraj-ul-Islam, and Abdul QM Khaliq. Local RBF method for multi- dimensional partial differential equations. Computers & Mathematics with Applica- tions, 74:292–324, 2017. [46] Siraj-ul-Islam and Imtiaz Ahmad. A comparative analysis of local meshless formu- lation for multi-asset option models. Engineering Analysis with Boundary Elements, 65:159–176, 2016. [47] WM Abd-Elhameed and YH Youssri. Spectral solutions for fractional differential equations via a novel Lucas operational matrix of fractional derivatives. Rom. J. Phys, 61(5-6):795–813, 2016. I. Ahmad et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5684 24 of 24 [48] Ihteram Ali, Sirajul Haq, Kottakkaran Sooppy Nisar, and Shams Ul Arifeen. Nu- merical study of 1D and 2D advection-diffusion-reaction equations using Lucas and Fibonacci polynomials. Arabian Journal of Mathematics, 10(3):513–526, 2021. [49] G Garmanjani, Roberto Cavoretto, and Mohsen Esmaeilbeigi. A RBF partition of unity collocation method based on finite difference for initial–boundary value prob- lems. Computers & Mathematics with Applications, 75(11):4066–4090, 2018. [50] Manzoor Hussain, Sirajul Haq, and Abdul Ghafoor. Meshless RBFs method for nu- merical solutions of two-dimensional high order fractional Sobolev equations. Com- puters & Mathematics with Applications, 79(3):802–816, 2020. [51] Peter D Lax. Weak solutions of nonlinear hyperbolic equations and their numerical computation. Communications on pure and applied mathematics, 7(1):159–193, 1954.