EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6838 ISSN 1307-5543 – ejpam.com Published by New York Business Global Tau Approach for the Time-Fractional Diffusion Equation Using Certain Chebyshev Polynomials W. M. Abd-Elhameed1,2, H. M. Alshehri2, M. H. Alharbi2, M. Adel3,∗, A. G. Atta4 1 Department of Mathematics, Faculty of Science, Cairo University, Giza 12613, Egypt 2 Department of Mathematics and Statistics, College of Science, University of Jeddah, Jeddah 21589, Saudi Arabia 3 Department of Mathematics, Faculty of Science, Islamic University of Madinah, Madinah, 42351, Saudi Arabia 4 Department of Mathematics, Faculty of Education, Ain Shams University, Cairo 11877, Egypt Abstract. This research article presents a numerical technique for solving the time-fractional diffusion equation (TFDE). The approach uses certain Chebyshev polynomials as basis functions. These polynomials are particular cases of the generalized Gegenbauer polynomials. The operational matrices of integer and fractional derivatives are used, along with the tau method, for the spatial and temporal discretization. Hence, the problem with its underlying conditions is converted into a system of equations that can be handled. We investigate the convergence of the double Chebyshev expansion and derive rigorous error bounds. Numerical examples show the method’s superior accuracy and efficiency over some existing methods in the literature. 2020 Mathematics Subject Classifications: 35R11, 65N35, 40A05 Key Words and Phrases: Time-fractional diffusion equation, Chebyshev polynomials, tau method, convergence analysis 1. Introduction The fractional differentiation extends the standard ordinary differentiation. Fractional differential equations (FDEs), in contrast to ordinary differential equations (DEs), are systems that include memory effects and long-range temporal dependencies that are well- represented by FDEs. This is particularly useful when they have to explain complex real-life occurrences. For instance, they model anomalous diffusion and stress-strain cor- relations in viscoelastic materials in the physical sciences. Using FDEs, we can better ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6838 Email addresses: waleed@cu.edu.eg (W. M. Abd-Elhameed), halshehri0271.stu@uj.edu.sa (H. M. Alshehri), mhhalharbi1@uj.edu.sa (M. H. Alharbi), adel@sci.cu.edu.eg (M. Adel), ahmed gamal@edu.asu.edu.eg (A. G. Atta) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 2 of 25 understand how past conditions impact current processes in biological systems, such as the spread of illness or changes in population size. Control theory, signal analysis, and electrical circuit modeling with memory components are other engineering fields that ben- efit from FDEs. They are also employed in environmental studies to represent the flow in a porous or fractured material and in finance to illustrate market linkages spanning several years. Due to their increased modeling capabilities and broad applicability, FDEs have been a hot topic in scientific computing and applied mathematics. Some applications of FDEs can be found in [1–4]. Due to the non-availability of analytic solutions for most FDEs, it is necessary to design effective numerical solutions for these equations. There are several numerical approaches for solving FDEs. The authors in [5] analyzed a numerical algorithm to treat the frac- tional Lane-Emden equations based on using the Laplace transform. The authors in [6] treated a nonlinear system of FDES utilizing the Adomian decomposition method. The authors in [7] used the differential transform method for analyzing fractional dynamical systems. In [8], the authors designed a numerical scheme to treat the fractional Burg- ers’ equation, which used the convolved Fermat polynomials as basis functions, enabling precise and computationally efficient solutions. The authors in [9] implemented a numeri- cal approach for the time-fractional generalized Kawahara equation using the eighth-kind CPs. In [10], a shifted Jacobi collocation scheme was developed to solve multidimensional time-fractional telegraph equations. In [11], the authors investigated a (3+1)-dimensional fractional coupled Burgers’ equation using a four-dimensional natural transform combined with the Adomian decomposition method. In addition, in [12], the authors presented a numerical scheme with adaptive step sizing to solve general FDEs. In [13], some analytical solutions of certain time-fractional heat order was treated. Modeling several physical phenomena across several fields still depends much on the time diffusion equation. Recent studies have used sophisticated computational methods and investigated new areas, therefore broadening their relevance. In computational math- ematics, several numerical approaches have been developed to solve the various forms of the diffusion equations more effectively. The authors in [14] proposed a kernel-based ap- proach to treat the TFDE. In [15], a finite difference–collocation scheme was designed for the TFDE, which combines the robustness of finite difference schemes with the advan- tages of using the collocation method. The authors in [16] used iterative methods for the TFDE. The authors in [17] offered a difference scheme for solving the TFDE. In [18], the authors analyzed a numerical approach to the TFDE. Hosoya polynomial was utilized in [19] to handle a class of TFDE, contributing an algebraic perspective to the numerical analysis. A numerical scheme based on B-splines was presented in [20] for addressing the TFDE. The authors in [21] used an efficient method for distributed-order TFDE. In [22], an exponential-sum-approximation technique was followed to tackle variable-order TFDE. The authors in [23] proposed analytical solutions for both time-fractional diffusion and convection–diffusion equations. Other contributions regarding some types of the TFDE can be found in [24–27]. Chebyshev polynomials (CPs) have vital roles in several branches of the applied sci- ences. In particular, they are extremely significant in areas like approximation theory, W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 3 of 25 numerical analysis, and applied mathematics. This versatility is due to their exceptional approximation properties; see, for example, [28, 29]. These polynomials are widely used in spectral and pseudospectral methods to solve various differential and integral equations. The most used kinds of CPs are the first and second kinds; see, for example, [30–33]. Other publications were devoted to the use of the third and fourth kinds. These four kinds are all particular ones of the classical Jacobi polynomials [34]. In addition to the classical Jacobi polynomials, there are other types of CPs. Masjed-Jamei, in his Ph.D. thesis [35], introduced two other trigonometric polynomials; he called them CPs of the fifth and sixth kinds. Various papers employed these polynomials; see, for example, [36, 37]. In addition, there are two other types of CPs, namely, CPs of the seventh and eighth kinds. They are specific polynomials of the generalized Gegenbauer polynomials. In [38], the seventh kind of CPs was used to solve the fractional diffusion equation, while the authors of [39] used the eighth kind of CPs to solve the nonlinear time-fractional generalized Kawahara equation. Other CPs were proposed and used to handle the FitzHugh-Nagumo equation in [40]. Spectral methods are numerical techniques employed to solve different types of differ- ential equations (DEs). The philosophy of applying this method is based on expanding the solution in terms of global basis functions such as trigonometric functions (Fourier series), orthogonal polynomials, or any other complete function systems. Spectral methods are able to attain exponential or “spectral” convergence rates, in contrast to the algebraic convergence that is produced by classic finite difference or finite element approaches that rely on local approximations. Because of this, they arise in fields like quantum physics, fluid dynamics, and wave propagation. For some applications of spectral methods, one can refer to [41–44] . Galerkin, collocation (or pseudospectral), and tau are the main spectral methods. The collocation approach, on the other hand, uses discrete places—the collocation points—to ensure that the differential equation is fulfilled. These collocation points are often chosen to represent the roots of orthogonal polynomials; see [45–49]. In the Galerkin method, we use basis functions and enforce the residual of the equation to be orthogonal to them; see, for example [50–54], The tau approach is more flexible than the Galerkin approach in choosing the basis functions, since there are no restrictions for the choice of them, and the underlying conditions are assumed to be as constraints; see, for example, [55–58]. The main objective of this article is to use certain CPs, which are particular polyno- mials of the generalized Gegenbauer polynomials, to solve the TFDE by applying the tau method. Theoretical and numerical studies of the errors obtained are given. The main contribution of the paper can be summarized in the following points: • Constructing the integer and fractional derivatives of the Chebyshev polynomials. • Applying the spectral tau method for spatial and temporal discretization. • Investigating the convergence and error analysis of the proposed expansion. • Testing the numerical algorithm by presenting some numerical results supported by some comparisons. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 4 of 25 The advantages of the proposed algorithm can be listed as follows: • Obtaining highly accurate approximations using a few terms of the basis functions. • Reducing the equation governed by its conditions into a system of equations that can be efficiently handled. • The capability of extending the tau approach for treating other types of FDEs. The paper follows the following structure: In the next section, we present the necessary mathematical preliminaries, including some properties of Caputo’s fractional derivative, and also some formulas concerned with the specific CPs that we will employ. Section 3 is interested in analyzing in detail the spectral Tau method for the TFDE using the proposed CPs. In Section 4, a thorough error analysis of the proposed Chebyshev expansion is con- ducted. Furthermore, we derive theoretical bounds for the double expansion. Section 5 ensures the applicability and accuracy of the presented algorithm through displaying sev- eral numerical experiments supported by some comparisons. Finally, concluding remarks are reported in Section 6. 2. Some fundamentals and formulas This section is confined to displaying an overview of the fractional calculus. Further- more, an account of the generalized Gegenbauer polynomials is given. Particular CPs that will be used as basis functions are introduced. Some of their formulas that are pivotal to deriving our proposed numerical algorithm are also provided. 2.1. Caputo’s fractional derivative Definition 1. [1], A fractional-order derivative, according to Caputo, is defined as Dν t ψ(s) = 1 Γ(k − ν) ∫ s 0 (s− t)k−ν−1ψ(k)(t)dt, ν > 0, s > 0, k − 1 < ν ≤ k, k ∈ N. (1) For Dν t with k − 1 < ν ≤ k, k ∈ N, the following identities hold: Dν t C =0, C is a constant, (2) Dν t s k = { 0, if k ∈ N0 and k < ⌈ν⌉, (k)! Γ(k−ν+1) s k−ν , if k ∈ N0 and k ≥ ⌈ν⌉, (3) where N = {1, 2, ...} and N0 = {0, 1, 2, . . .}, and ⌈ν⌉ is the ceiling function. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 5 of 25 2.2. An account on generalized Gegenbauer polynomials We give an account of the generalized Gegenbauer polynomials G (µ,β) k (x). They are orthogonal on [−1, 1] regarding w(x) = (1 − x2)µ− 1 2 |x|2β. They can be represented as [59, 60] G (µ,β) k (x) =  (µ+ β) k 2( β + 1 2 ) k 2 P (µ− 1 2 ,β− 1 2) k 2 ( 2x2 − 1 ) , if k is even, (µ+ β) k+1 2( β + 1 2 ) k+1 2 xP (µ− 1 2 ,β+ 1 2) k−1 2 ( 2x2 − 1 ) , if k is odd, (4) where (z)k = Γ(z+k) Γ(z) denotes the Pochhammer symbol and P (µ,β) k (x) are the standard Jacobi polynomials. The orthogonality relation for G (µ,β) k (x) is 1∫ −1 w(x)G (µ,β) k (x)G(µ,β) m (x) dx = { hµ,βk , k = m, 0, k ̸= m, (5) where hµ,βk is given by hµ,βk = ( Γ ( β + 1 2 ))2 (k + µ+ β) (Γ(µ+ β))2  Γ ( µ+ k+1 2 ) Γ ( k 2 + µ+ β )( k 2 ) ! Γ ( 1+k 2 + β ) , k even, Γ ( k 2 + µ ) Γ ( k 2 + µ+ β + 1 2 )( k−1 2 ) ! Γ ( k 2 + β + 1 ) , k odd. (6) Remark 1. Several particular polynomials, including CPs, are particular ones of the generalized polynomials G (µ,β) k (x), see [40]. The authors in the same paper introduced certain CPs as particular ones of G (µ,β) k (x), and derived some formulas regarding them. 2.3. Some properties and formulas of a certain kind of CPs The authors in [40] proposed the polynomials Nk(x) = G2,1 k (x). They are defined as Nk(x) =  (3) k 2 ( 3 2) k 2 P ( 3 2 , 1 2) k 2 ( 2x2 − 1 ) , if k even, x(3) k+1 2 ( 3 2) k+1 2 P ( 3 2 , 3 2) k−1 2 ( 2x2 − 1 ) , if k odd. (7) {Nk(x)}k≥0 are orthogonal on [−1, 1] in the sense that∫ 1 −1 x2 ( 1− x2 )3/2 Nk(x) Nj(x) dx = hk, (8) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 6 of 25 where hn = π 128  (k + 2) (k + 4), if k = j, k even, (k + 1) (k + 5), if k = j, k odd, 0, if k ̸= j. (9) The following two lemmas give the explicit analytic form ofNk(x) and its inversion formula. Lemma 1. [40] For every positive integer s, the polynomials Ns(x) have the following expression: Ns(x) = ⌊ s 2⌋∑ m=0 (−1)m (3)s−m m! (⌊ s 2 ⌋ −m ) ! ( 3 2 ) −m+⌊ 1+s 2 ⌋ xs−2m. (10) Lemma 2. [40] For every positive integer s, xs can be expanded as xs = ⌊ s 2⌋∑ m=0 (s− 2m+ 3) ⌊ s 2 ⌋ ! ( 3 2 ) ⌊ s+1 2 ⌋ m! (3)s−m+1 Ns−2m(x). (11) Now, we give the high-order derivatives expression for the polynomials G (µ,β) k (x) in terms of their original polynomials. Theorem 1. [40] The m-th derivative: dk Ns(x) d xk can be represented as dmNs(x) d xm = ⌊ s−m 2 ⌋∑ r=0 Gm r,sNs−m−2r(x), s ≥ m, (12) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 7 of 25 where Gm r,s = 2m  (1 + s−m)(3 + s− 2r −m)(s+ 2)! (s+ 1)!r!(s− r −m+ 3)! × 4F3 ( −r,−1 2 − s 2 , 1 2 − s 2 + m 2 ,−3− s+ r +m −2− s, 12 − s 2 ,− 1 2 − s 2 + m 2 ∣∣∣∣∣ 1 ) , s even,m even, (−1− s+m)(−3− s+ 2r +m)(s+ 1)! r! (s− r −m+ 3)! × 4F3 ( −r,−1− s 2 , 1 2 − s 2 + m 2 ,−3− s+ r +m −2− s,− s 2 ,− 1 2 − s 2 + m 2 ∣∣∣∣∣ 1 ) , s odd,m odd, (−2− s+m)(−3− s+ 2r +m)(s+ 1)! r! (s− r −m+ 3)! × 4F3 ( −r,−1− s 2 ,− s 2 + m 2 ,−3− s+ r +m −2− s,− s 2 ,−1− s 2 + m 2 ∣∣∣∣∣ 1 ) , s odd,m even, (2 + s−m)(3 + s− 2r −m)(s+ 2)! (1 + s) r! (s− r −m+ 3)! × 4F3 ( −r,−1 2 − s 2 ,− s 2 + m 2 ,−3− s+ r +m −2− s, 12 − s 2 ,−1− s 2 + m 2 ∣∣∣∣∣ 1 ) , s even,m odd. (13) 2.4. The shifted CPs Now, we introduce the corresponding shifted polynomials on [0, 1] to the polynomials Nk(x). We denote them by H∗ k(x) and define them as H∗ k(x) = Nk(2x− 1). (14) All the formulas of the polynomials Nk(x) can be transformed to give their counterparts on [0, 1]. From (8), the set {H∗ k(x)}k≥0 is an orthogonal set on [0, 1] with∫ 1 0 H∗ n(x)H∗ m(x)w∗(x) dx = 1 16 hn, (15) with w∗(x) = (1− 2x)2(x(1− x))3/2, and hn is defined in (9). Theorem 2. The mth-derivative of H∗ j (x) can be expressed as dmH∗ s(x) d xm = ⌊ s−m 2 ⌋∑ r=0 B̃m r,sH∗ s−m−2r(x), s ≥ m, (16) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 8 of 25 with B̃k r,j = 2k Gk r,j , where Gk r,j is given in (13). Proof. Replacing x by (2x− 1) in Theorem 1, we get the result in (16). 3. Tau procedure for the TFDE Consider the following TFDE [20, 61]: Dζ t λ(x, t)− β λxx(x, t) = g(x, t), 0 < ζ ≤ 1, (17) governed by the following conditions: λ(x, 0) = ρ1(x), 0 < x < 1, (18) λ(0, t) = ρ2(t), λ(1, t) = ρ3(t), 0 < t < τ, (19) where β is a positive constant, ρ1(x), ρ2(t), and ρ3(t) are known continuous functions, and g(x, t) is the source term. If we let PM = span{H∗ i (x)H∗ j (t) : 0 ≤ i, j ≤ M}, then; any function λM(x, t) ∈ PM can be represented as λM(x, t) = M∑ i=0 M∑ j=0 cij H∗ i (x)H∗ j (t) = H∗(x)CH∗(t)T , (20) where H∗(t)T = [H∗ 0(t),H∗ 1(t), . . . ,H∗ M(t)]T , and C = (cij)0≤i,j≤M is the (M + 1)2- dimensional matrix of unknowns. The residual RM(x, t) of Eq (17) is given by RM(x, t) = Dζ t λ M(x, t)− β λMxx(x, t)− g(x, t). (21) If we apply the tau method, then we get∫ 1 0 ∫ 1 0 RM(x, t)H∗ r(x)H∗ s(t))w ∗(x)w∗(t) dx dt = 0, 1 ≤ r ≤ M− 1, 1 ≤ s ≤ M. (22) Now, consider the following matrices: G = (gr,s)(M−1)×M, gr,s = ∫ 1 0 ∫ 1 0 g(x, t)H∗ r(x)H∗ s(t)w ∗(x)w∗(t) dx dt, (23) Z = (zi,r)(M+1)×(M−1), zi,r = ∫ 1 0 H∗ i (x)H∗ r(x)w ∗(x) dx, (24) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 9 of 25 B = (bj,s)(M+1)×M, bj,s = ∫ 1 0 H∗ j (t)H∗ s(t)w ∗(t) dt, (25) Y = (yi,r)(M+1)×(M−1), yi,r = ∫ 1 0 d2H∗ i (x) d x2 H∗ r(x)w ∗(x) dx, (26) K = (kj,s)(M+1)×M, kj,s = ∫ 1 0 Dζ tH∗ j (t)H∗ s(t)w ∗(t) dt. (27) Therefore, Eq. (22) can be rewritten as M∑ i=0 M∑ j=0 cij zi,r kj,s − β M∑ i=0 M∑ j=0 cij yi,r bj,s = fr,s, 0 ≤ r ≤ M− 2, 0 ≤ s ≤ M− 1, (28) or in the following matrix form: ZT CK− βYTC B = G. In addition, the conditions in (18) and (19) lead to the following equations: M∑ i=0 M∑ j=0 cij ai,r H∗ j (0) = ∫ 1 0 ρ1(x)H∗ r(x)w ∗(x) dx, 0 ≤ r ≤ M, (29) M∑ i=0 M∑ j=0 cij aj,sH∗ i (0) = ∫ 1 0 ρ2(t)H∗ s(t)w ∗(t) dt, 0 ≤ s ≤ M− 1, (30) M∑ i=0 M∑ j=0 cij aj,sH∗ i (1) = ∫ 1 0 ρ3(t)H∗ s(t)w ∗(t) dt, 0 ≤ s ≤ M− 1. (31) The resultant algebraic system of equations, consisting of Eqs. (28)-(31), has an order of (M+ 1)2 and may be solved using an appropriate method. Now, we give an explicit form for the elements of the matrices Z,K,Y . Theorem 3. The elements zi,r, bj,s, yi,r, and kj,s have the following forms: zi,r = ∫ 1 0 H∗ i (x)H∗ r(x)w ∗(x) dx = 1 16 hi, (32) bj,s = ∫ 1 0 H∗ j (t)H∗ s(t)w ∗(t) dt = 1 16 hj , (33) yi,j = ∫ 1 0 d2H∗ i (x) dθ2 H∗ j (x)w ∗(x) dx = 1 16 ⌊ i−2 2 ⌋∑ r=0 B̃2 r,i hi−k−2r, (34) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 10 of 25 and kj,s = ∫ 1 0 Dζ tH∗ j (t)H∗ s(t)w ∗(t) dt = j∑ r=1 j∑ m=0 s∑ p=0 s∑ n=0 2rr!(−1)) j−m 2 (−1)m−ra(j +m) ( m r ) (3) j+m 2( j−m 2 ) !Γ(r − α+ 1) ( 3 2 ) m−j 2 +⌊ j+1 2 ⌋ Γ ( 1 2(−j +m+ 2) + ⌊ j 2 ⌋)× 3 √ π2p−2(−1)n−p(−1) s−n 2 a(n+ s) ( n p ) (3)n+s 2 ( α2 − α(2p+ 2r + 1) + (p+ r)2 + p+ r + 5 ) s−n 2 ! ( 3 2 ) n−s 2 +⌊ s+1 2 ⌋ Γ ( 1 2(n− s+ 2) + ⌊ s 2 ⌋) × Γ ( p+ r − α+ 5 2 ) Γ(p+ r − α+ 7) . (35) Proof. With the direct application of the orthogonality relation (15), we can easily acquire the elements zi,r and bj,s. With the direct application of (15) with (16), we can easily obtain the elements yi,j . Now, we will compute kj,s, Based on Eq. (10), the power form representation of Nj(t) as can be rewritten as Nj(t) = j∑ r=0 (−1) j−r 2 aj+r (3)j− j−r 2( j−r 2 ) ! ( 3 2 ) ⌊ j+1 2 ⌋− j−r 2 Γ ( −1 2(j − r) + ⌊ j 2 ⌋ + 1 ) tr, (36) where ar = { 1, if r even, 0, otherwise. (37) and hence, H∗ j (t) given in (14) may be represented as H∗ j (t) = j∑ r=0 (−1) j−r 2 aj+r (3)j− j−r 2( j−r 2 ) ! ( 3 2 ) ⌊ j+1 2 ⌋− j−r 2 Γ ( −1 2(j − r) + ⌊ j 2 ⌋ + 1 ) (2 t− 1)r, (38) the last equation can be rewritten after using the relation (2 t−1)2 = r∑ n=0 (2 t)n(−1)r−n ( r n ) , expanding, rearranging and collecting similar terms as H∗ j (t) = j∑ r=0 j∑ m=0 2r (−1) j−m 2 (−1)m−r aj+m ( m r ) (3) j+m 2( j−m 2 ) ! ( 3 2 ) m−j 2 +⌊ j+1 2 ⌋ Γ ( 1 2(−j +m+ 2) + ⌊ j 2 ⌋) tr. (39) Therefore, the identity of fractional Caputo derivative (3) along with Eq. (39) enables us W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 11 of 25 to write Dζ tH∗ j (t) = j∑ r=1 j∑ m=0 2rr!(−1) j−m 2 (−1)m−r aj+m ( m r ) (3) j+m 2( j−m 2 ) ! Γ(r − α+ 1) ( 3 2 ) m−j 2 +⌊ j+1 2 ⌋ Γ ( 1 2(−j +m+ 2) + ⌊ j 2 ⌋) tr−ζ . (40) Now, kj,s is given by kj,s = ∫ 1 0 Dζ tH∗ j (t)H∗ s(t)w ∗(t) dt = j∑ r=1 j∑ m=0 2rr!(−1) j−m 2 (−1)m−r aj+m ( m r ) (3) j+m 2( j−m 2 ) ! Γ(r − α+ 1) ( 3 2 ) m−j 2 +⌊ j+1 2 ⌋ Γ ( 1 2(−j +m+ 2) + ⌊ j 2 ⌋)× ∫ 1 0 tr−ζH∗ s(t)w ∗(t) dt = j∑ r=1 j∑ m=0 2rr!(−1) j−m 2 (−1)m−r aj+m ( m r ) (3) j+m 2( j−m 2 ) ! Γ(r − α+ 1) ( 3 2 ) m−j 2 +⌊ j+1 2 ⌋ Γ ( 1 2(−j +m+ 2) + ⌊ j 2 ⌋) × s∑ p=0 s∑ n=0 2p(−1)n−pis−na(n+ s) ( n p ) (3)n+s 2 s−n 2 ! ( 3 2 ) n−s 2 +⌊ s+1 2 ⌋ Γ ( 1 2(n− s+ 2) + ⌊ s 2 ⌋) ∫ 1 0 (1− 2 t)2(t (1− t))3/2t−α+p+r dt. (41) Finally, after evaluating the integral on the right-hand side of the previous equation, we get the desired result (35). 4. Error bound In this part, we analyze the errors in the proposed polynomial expansion in detail. It is planned to prove four theorems. • The first theorem provides an upper bound for the truncation error. • The second theorem, provides a maximum estimate for the m-th derivative of trun- cation error with respect to the variable x. • The third theorem provides a maximum estimate for the fractional derivative of the truncation error with respect to the variable t • The last theorem shows ∥RM(x, t)∥L2 ω(Ω) for sufficiently high M, will be small enough, where ω = w∗(x)w∗(t). Lemma 3. [62] Given that m ≥ 1 and that m+c > 1 and m+d > 1, respectively, for all constants c,d, one has Γ(m + c) Γ(m + d) ≤ oc,dm mc−d, (42) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 12 of 25 where oc,dm = exp ( c− d 2 (m + d− 1) + 1 12 (m + c− 1) + (c− d)2 m ) . (43) Remark 2. For fixed c,d. oc,dm can be written as: oc,dm = 1 +O(m−1). Theorem 4. Assume that ∂i+j λ(x,t) ∂ xi ∂ tj ∈ C(Ω), i, j = 0, 1, 2, . . . ,M+ 1 and λM(x, t) is the suggested approximate solution belonging to PM and ℓM = sup (x,t)∈Ω ∣∣∣∣∣∂2 (M+1) λ(x, t) ∂ xM+1 ∂ tM+1 ∣∣∣∣∣ . (44) Consequently, this estimate is valid: ∥λ(x, t)− λM(x, t)∥L2 ω(Ω) ≲ ℓM M 1 2 ((M+ 1)!)2 , (45) where q1 ≲ q2 implies the existence of a constant n satisfies q1 ≤ n q2. Proof. Consider the following Taylor expansion of λ(x, t): χM(x, t) = M∑ i=0 M−i∑ j=0 ( ∂i+j λ(x, t) ∂ xi ∂ tj ) (0,0) xi tj i! j! , (46) λ(x, t)− χM(x, t) = xM+1 tM+1 ∂2 (M+1) λ(n1, n2) ((M+ 1)!)2 ∂ xM+1 ∂ tM+1 , (n1, n2) ∈ Ω. (47) Since zM(x, t) is the best approximate solution of z(x, t), we get, in accordance with the best approximation concept: ∥λ(x, t)− λM(x, t)∥2L2 ω(Ω) ≤ ∥λ(x, t)− χM(x, t)∥2L2 ω(Ω) = ∫ 1 0 ∫ 1 0 ℓM 2 x2 (M+1) t2 (M+1) ((M+ 1)!)4 ω dx d t = ℓM 2 π ((M+ 1)!)4 ( Γ ( 2M+ 9 2 ) Γ(2M+ 5) − 5 Γ ( 2M+ 11 2 ) Γ(2M+ 6) + 8 Γ ( 2M+ 13 2 ) Γ(2M+ 7) −4 Γ ( 2M+ 15 2 ) Γ(2M+ 8) )2 . (48) The following estimate may be written using Lemma 3. ∥λ(x, t)− λM(x, t)∥L2 ω(Ω) ≲ ℓM M 1 2 ((M+ 1)!)2 . (49) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 13 of 25 Theorem 5. Suppose that λ(x, t), λM(x, t) and ∂i+j λ(x,t) ∂ xi ∂ tj meet the assumption of Theorem 4 and τM,m = sup (x,t)∈Ω ∣∣∣∣ ∂2M−m+2 λ(x, t) ∂ xM−m+1 ∂ tM+1 ∣∣∣∣ , m ∈ N, (50) where N = {1, 2, . . .}. Then, we can write∥∥∥∥ ∂m ∂ xm (λ(x, t)− λM(x, t)) ∥∥∥∥ L2 ω(Ω) ≲ τM,m (M (M−m)) 1 4 (M+ 1)! (M−m+ 1)! . (51) Proof. Assume that ∂m χM(x,t) ∂ xm is the Taylor expansion of ∂m λ(x,t) ∂ xm about the point (0, 0), then the residual between ∂m λ(x,t) ∂ xm and ∂m χM(x,t) ∂ xm can be written as [63] ∂m ∂ xm (λ(x, t)− χM(x, t)) = xM−m+1 tM+1 ∂2M−m+2 λ(n̄1, n̄2) (M+ 1)! (M−m+ 1)! ∂ xM−m+1 ∂ tM+1 , (n̄1, n̄2) ∈ Ω. (52) Since ∂m λM(x,t) ∂ xm is the best approximate solution of ∂m λ(x,t) ∂ xm , then according to the defini- tion of the best approximation, we get∥∥∥∥ ∂m ∂ xm (λ(x, t)− λM(x, t)) ∥∥∥∥ L2 ω(Ω) ≤ ∥∥∥∥ ∂m ∂ xm (λ(x, t)− χM(x, t)) ∥∥∥∥ L2 ω(Ω) . (53) We get the desired result by following similar steps to those followed in Theorem 4. Theorem 6. Suppose that Dζ t λ(x, t) ∈ C(Ω) satisfies the conditions of Theorem 4, then the following estimation holds∥∥∥Dζ t (λ(x, t)− λM(x, t)) ∥∥∥ L2 ω(Ω) ≲ ℓM M 1 4 (M− ζ) 1 4 (M+ 1)! Γ(M+ 2− ζ) . (54) Proof. According to Eq. (47) and properties of the Caputo operator in (3), one gets∣∣∣Dζ t (λ(x, t)− χM(x, t) ) ∣∣∣ ≤ xM+1 tM+1−ζ ℓM Γ(M+ 2− ζ) (M+ 1)! , (n1, n2) ∈ Ω. (55) Taking ∥.∥L2 ω(Ω) yields∥∥∥Dζ t (λ(x, t)− λM(x, t)) ∥∥∥ L2 ω(Ω) ≤ ∫ 1 0 ∫ 1 0 ℓM 2 x2 (M+1) t2 (M+1−ζ) (Γ(M+ 2− ζ))2((M+ 1)!)2 ω dx d t. (56) Now, imitating similar steps as in Theorem 4, we get∥∥∥Dζ t (λ(x, t)− λM(x, t)) ∥∥∥ L2 ω(Ω) ≲ ℓM M 1 4 (M− ζ) 1 4 (M+ 1)! Γ(M+ 2− ζ) . (57) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 14 of 25 Theorem 7. Assume that RM(x, t) be the residual of Eq. (17), then ∥RM(x, t)∥L2 ω(Ω) will be sufficiently small for the sufficiently large values of M. Proof. Eqs. (21) and (17) enables us too write RM(x, t) as RM(x, t) = Dζ t ( λM(x, t)− λ(x, t) ) − β ∂2 ∂ x2 ( λM(x, t)− λ(x, t) ) . (58) If we consider L2-norm, then using Theorems 4,5 and 6, we get ∥RM(x, t)∥L2 ω(Ω) ≲ ℓM M 1 4 (M− ζ) 1 4 (M+ 1)! Γ(M+ 2− ζ) − β τM,2 (M (M− 2)) 1 4 (M+ 1)! (M− 2 + 1)! . (59) For large enough values of M, it is evident from Eq. (59) that ∥RM(x, t)∥L2 ω(Ω) will be small enough. This concludes the proof of the theorem. 5. Some numerical tests In this section, we present some numerical examples to demonstrate the accuracy and efficiency of the suggested numerical algorithm by using the absolute error (AE), maximum absolute error (MAE) and L∞ - error AE = ∣∣λ(x, t)− λM(x, t) ∣∣ , (60) MAE = max xi ∣∣λ(xi, xj)− λM(xi, xj) ∣∣ , xi ∈ { 1 10 , 2 10 , 3 10 , . . . , 1 } , (61) L∞ = max (x,t)∈Ω ∣∣λ(x, t)− λM(x, t) ∣∣ . (62) In addition, comparisons with some methods are given. Example 1. [64] Consider the following equation: Dζ t λ(x, t)− λxx(x, t) = 1 Γ(2− α) t1−α x2 (1− x) + 2 t (3x− 1), 0 < α ≤ 1, (63) subject to the initial condition (IC) λ(x, 0) = 0, 0 < x ≤ 1, (64) and the homogeneous boundary conditions (HBCs) λ(0, t) = λ(1, t) = 0, 0 < t ≤ 1, (65) where λ(x, t) = t x2 (1− x) is the exact solution of this problem. Table 1 shows the AE at various values of ζ when M = 3. Figure 1 illustrates the AE W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 15 of 25 Table 1: The AE of Example 1 at M = 3. (x, t) ζ = 0.2 CPU time ζ = 0.4 CPU time ζ = 0.6 CPU time ζ = 0.8 CPU time (0.1,0.1) 3.74914× 10−14 8.69194× 10−15 8.28406× 10−15 1.75789× 10−14 (0.2,0.2) 5.27894× 10−14 1.66768× 10−14 1.26583× 10−14 2.83619× 10−14 (0.3,0.3) 6.04031× 10−15 9.24955× 10−15 7.10196× 10−15 1.20529× 10−14 (0.4,0.4) 6.2942× 10−14 1.40582× 10−14 6.8695× 10−16 1.98938× 10−14 (0.5,0.5) 1.13222× 10−13 2.469 4.27158× 10−14 2.656 1.77636× 10−15 2.999 4.69347× 10−14 2.578 (0.6,0.6) 1.07289× 10−13 6.05765× 10−14 9.96425× 10−15 5.43871× 10−14 (0.7,0.7) 5.93969× 10−14 5.37348× 10−14 1.57374× 10−14 4.09117× 10−14 (0.8,0.8) 9.90874× 10−15 2.01922× 10−14 8.20177× 10−15 1.8055× 10−14 (0.9,0.9) 3.08087× 10−15 1.78052× 10−14 1.01724× 10−14 2.45637× 10−15 at ζ = 0.5 and ζ = 0.9 when M = 3. Also, the CPU time (in seconds) for the method is computed in this table. Table 2 gives a comparison of L2 and L∞ errors between our method at M = 3 and the method in [64] at ∆t = 0.001 and ∆x = 0.015625. Figure 1: The AE for Example 1 at ζ = 0.5 and ζ = 0.9 when M = 3. Table 2: Comparison of errors for Problem 1. Method in [64] at ∆t = 0.001 and ∆x = 0.015625 Our method at M = 3 L2 error L∞ error L2 error L∞ error 1.6868× 10−12 2.3978× 10−12 6.37117× 10−16 1.44382× 10−13 Example 2. [20] Consider the following equation: Dζ t λ(x, t)− λxx(x, t) = g(x, t), 0 < ζ ≤ 1, (66) controlled by λ(x, 0) = 0, 0 < x < 1, λ(0, t) = λ(1, t) = 0, 0 < t < 1, (67) where g(x, t) is chosen to meet the exact solution given by λ(x, t) = sin(π t) sin(π x). Table 3 gives a comparison of MAE between our method at M = 8 and the method in [20] W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 16 of 25 at M = 32 and ∆x = 0.001. Figures 2 and 3 illustrate the MAE and L∞-error at various values of M when ζ = 0.5. Table 4 shows the AE at various values of t when ζ = 0.7 and M = 8. Table 5 shows the MAE and L∞-error at various values of M when ζ = 0.9. Also, Table 6 shows the MAE and L∞-error at various values of M when ζ = 0.3. Table 3: Comparison of the MAE for Example 2. ζ Method in [20] (M = 32 and ∆x = 0.001) Our method (M = 8) 0.5 1.10× 10−3 3.30215× 10−6 0.7 3.21× 10−3 3.30734× 10−6 ◆ ◆ ◆ ◆ ◆ ◆ ◆ 2 3 4 5 6 7 8 10 5 10 4 0.001 0.010 0.100 ℳ E rr o rs ◆ MAEs Figure 2: The MAE of Example 2 at various values of M when ζ = 0.5. Table 4: The AE of Example 2 at ζ = 0.7, M = 8. x t = 0.2 t = 0.4 t = 0.6 t = 0.8 0.1 2.98478× 10−7 5.82311× 10−7 6.6952× 10−7 4.94957× 10−7 0.2 6.22613× 10−7 1.19821× 10−6 1.36547× 10−6 9.99651× 10−7 0.3 9.18475× 10−7 1.75345× 10−6 1.98639× 10−6 1.44529× 10−6 0.4 1.33961× 10−6 2.48775× 10−6 2.76405× 10−6 1.96795× 10−6 0.5 1.57757× 10−6 2.89203× 10−6 3.18343× 10−6 2.24218× 10−6 0.6 1.33959× 10−6 2.488× 10−6 2.76412× 10−6 1.96804× 10−6 0.7 9.18451× 10−7 1.75373× 10−6 1.98646× 10−6 1.44539× 10−6 0.8 6.22605× 10−7 1.19833× 10−6 1.36551× 10−6 9.99693× 10−7 0.9 2.98475× 10−7 5.82369× 10−7 6.69545× 10−7 4.94976× 10−7 W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 17 of 25 2 3 4 5 6 7 8 10 -5 10 -4 0.001 0.010 0.100 ℳ E rr o rs L∞-error Figure 3: The L∞-error of Example 2 at various values of M when ζ = 0.5. Table 5: The errors of Example 2 at various values of M when ζ = 0.9. M 2 4 6 8 MAE 2.90874×10−1 1.55707×10−2 2.98468×10−4 3.31145×10−6 L∞-error 2.91952×10−1 1.57822×10−2 3.02187×10−4 3.35344×10−6 Table 6: The errors of Example 2 at various values of M when ζ = 0.3. M 2 3 4 5 6 7 8 MAE 3.0980×10−1 2.854×10−1 2.8600×10−2 1.4596×10−2 3.0076×10−4 2.8939×10−4 3.3021×10−6 L∞-error 3.0986×10−1 2.8600×10−1 1.5560×10−2 1.4612×10−2 3.0092×10−4 2.8949×10−4 3.3090×10−6 Example 3. Consider the following equation: Dζ t λ(x, t)− λxx(x, t) = ( π2 t2 + 2 t2−α Γ(3− α) ) sin(π x), 0 < α ≤ 1, (68) subject to the IC λ(x, 0) = 0, 0 < x ≤ 1, (69) and the HBCs λ(0, t) = λ(1, t) = 0, 0 < t ≤ 1, (70) where λ(x, t) = t2 sin(π x) is the exact solution of this problem. Figure 4 illustrates the MAE and L∞-error at various values of M when ζ = 0.5. Table 7 shows the AE at various values of t when ζ = 0.9 and M = 8. Table 8 shows the MAE and L∞-error at various values of M when ζ = 0.9. Figure 5 illustrates the AE at various values of M when ζ = 0.3. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 18 of 25 ◆ ◆ ◆ ◆ ◆ ◆ ◆ 2 3 4 5 6 7 8 10 6 10 5 10 4 0.001 0.010 0.100 ℳ E rr o rs ◆ MAEs L∞-error Figure 4: The errors of Example 3 at various values of M when ζ = 0.5. Table 7: The AE of Example 3 at ζ = 0.9, M = 8. x t = 0.2 t = 0.4 t = 0.6 t = 0.8 0.1 1.1953× 10−8 6.96178× 10−8 1.83032× 10−7 3.5209× 10−7 0.2 2.62526× 10−8 1.47194× 10−7 3.81809× 10−7 7.30076× 10−7 0.3 3.98826× 10−8 2.18594× 10−7 5.63519× 10−7 1.44529× 10−6 0.4 6.41293× 10−8 3.26744× 10−7 8.22244× 10−7 1.07306× 10−6 0.5 7.8752× 10−8 3.89168× 10−7 9.68425× 10−7 1.54566× 10−6 0.6 6.41292× 10−8 3.26811× 10−7 8.22324× 10−7 1.80976× 10−6 0.7 3.98859× 10−8 2.1865× 10−7 5.63611× 10−7 1.54565× 10−6 0.8 2.62589× 10−8 1.47195× 10−7 3.81857× 10−7 7.3006× 10−7 0.9 1.19538× 10−8 6.96403× 10−8 1.83059× 10−7 3.52091× 10−7 Table 8: The errors of Example 3 at various values of M when ζ = 0.7. M 2 4 6 8 MAE 8.7641×10−2 4.11803×10−3 7.9429×10−5 8.71518×10−7 L∞-error 2.54328×10−1 1.32×10−2 2.66423×10−4 2.94457×10−6 Example 4. [64] Consider the following equation Dζ t λ(x, t)− λxx(x, t) = ( 4π2 t2 + 2 t2−α Γ(3− α) ) sin(2π x), 0 < α ≤ 1, (71) subject to the Initial condition (IC): λ(x, 0) = 0, 0 < x ≤ 1, (72) W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 19 of 25 Figure 5: The AE of Example 3 at various values of M when ζ = 0.3. and the HBCs λ(0, t) = λ(1, t) = 0, 0 < t ≤ 1, (73) where λ(x, t) = t2 sin(2π x) is the exact solution of this problem. Table 9 gives a comparison of MAE between our method at M = 9 and method in [64] at ∆t = 0.001 and M = 64. Table 10 displays the AE at M = 9 and ζ = 0.5. Table 11 reports the MAE and L∞-error at various values of M when ζ = 0.8 Figure 6 illustrates the AE at various values of t when M = 9 and ζ = 0.3. Table 9: Comparison of the MAE for Example 4. Method in [64] at ∆t = 0.001 and M = 64 Our method at M = 9 7.70× 10−4 1.66× 10−5 Table 11: The errors of Example 4 at various values of M when ζ = 0.8. M 3 6 9 MAEs 1.37846×10−1 1.42851×10−2 1.5862×10−5 L∞-error 2.83552×10−1 3.13749×10−2 4.61017×10−5 W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 20 of 25 Table 10: The AE of Example 4 at M = 9 and ζ = 0.5. x t = 0.2 t = 0.4 t = 0.6 t = 0.8 0.1 4.33573× 10−7 1.71139× 10−6 4.02304× 10−6 7.16437× 10−6 0.2 7.58102× 10−7 3.01934× 10−6 7.05377× 10−6 1.25795× 10−5 0.3 1.16957× 10−6 4.53537× 10−6 1.06733× 10−5 1.89359× 10−5 0.4 1.86083× 10−6 7.1667× 10−6 1.66918× 10−5 2.95361× 10−5 0.5 2.46676× 10−8 9.86345× 10−8 5.44645× 10−8 1.4437× 10−8 0.6 1.81921× 10−6 7.3363× 10−6 1.66× 10−5 2.95617× 10−5 0.7 1.1434× 10−6 4.64781× 10−6 1.06156× 10−5 1.89542× 10−5 0.8 7.42571× 10−7 3.08815× 10−6 7.01915× 10−6 1.25912× 10−5 0.9 4.25343× 10−7 1.74765× 10−6 4.00476× 10−6 7.17045× 10−6 0.0 0.2 0.4 0.6 0.8 1.0 0.00000 0.00001 0.00002 0.00003 0.00004 x E rr o rs t=0.1 t=0.3 t=0.5 t=0.7 Figure 6: The AE of Example 4 at M = 9 and ζ = 0.3. 6. Concluding remarks This article presented an effective Tau-based numerical scheme for solving the TFDE based on employing certain CPs, which are particular polynomials of the generalized Gegenbauer polynomials. Our tau algorithm transforms the equation with its initial- boundary equations into an algebraic system of equations that can be numerically solved. The method exhibits excellent approximation capabilities, as evidenced by the low errors in various test problems. The approach used in this paper may be a powerful tool for handling other fractional PDEs in future research. In addition, we anticipate that other generalized CPs as special cases of the generalized Gegenbauer polynomials may be introduced and used along with suitable spectral methods to treat different types of differential equations. Acknowledgements The researchers wish to extend their sincere gratitude to the Deanship of Scientific Research at the Islamic University of Madinah for the support provided to the Post- Publishing Program. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 21 of 25 Declarations • Conflicts of interest: The authors declare that they have no conflicts of interest. References [1] I. Podlubny. Fractional Differential Equations: An Introduction to Fractional Deriva- tives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Elsevier, 1998. [2] R. Hilfer. Applications of Fractional Calculus in Physics. World Scientific, Singapore, 2000. [3] M.D. Ortigueira. Fractional Calculus in Applied Sciences and Engineering. Springer, Cham, 2015. [4] V.E. Tarasov. Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields, and Media. Springer, Berlin, 2011. [5] R. Saadeh, A. Burqan, and A. El-Ajou. Reliable solutions to fractional Lane- Emden equations via Laplace transform and residual error function. Alex. Eng. J., 61(12):10551–10562, 2022. [6] A. Afreen and A. Raheem. Study of a nonlinear system of fractional differential equations with deviated arguments via Adomian decomposition method. Int. J. Appl. Comput. Math., 8(5):269, 2022. [7] A. Rysak and M. Gregorczyk. Differential transform method as an effective tool for investigating fractional dynamical systems. Appl. Sci., 11(15):6955, 2021. [8] W.M. Abd-Elhameed, O.M. Alqubori, N.M.A. Alsafri, A.K. Amin, and A.G. Atta. A matrix approach by convolved Fermat polynomials for solving the fractional Burgers’ equation. Mathematics, 13(7):1135, 2025. [9] H.M. Ahmed, R.M. Hafez, and W.M. Abd-Elhameed. A computational strategy for nonlinear time-fractional generalized Kawahara equation using new eighth-kind Chebyshev operational matrices. Phys. Scr., 99(4):045250, 2024. [10] R.M. Hafez and Y.H. Youssri. Shifted Jacobi collocation scheme for multidimensional time-fractional order telegraph equation. Iran. J. Numer. Anal. Optim., 10(1):195– 223, 2020. [11] H. Alsaud and H. Eltayeb. The four-dimensional natural transform Adomian decom- position method and (3+1)-dimensional fractional coupled Burgers’ equation. Frac- tals, 8(4):227, 2024. [12] M. Shams and A. Alalyani. High performance adaptive step size fractional numerical scheme for solving fractional differential equations. Sci. Rep., 15(1):13006, 2025. [13] A.M.S. Mahdy, Kh. Lotfy, E.A. Ismail, A. El-Bary, M. Ahmed, and A.A. El-Dahdouh. Analytical solutions of time-fractional heat order for a magneto-photothermal semi- conductor medium with Thomson effects and initial stress. Results Phys., 18:103174, 2020. [14] M. Fardi. A kernel-based method for solving the time-fractional diffusion equation. Numer. Methods Partial Differ. Equations, 39(3):2719–2733, 2023. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 22 of 25 [15] S. Kumar, R.K. Pandey, K. Kumar, S. Kamal, and T.N. Dinh. Finite difference– collocation method for the generalized fractional diffusion equation. Fractal Fract., 6(7):387, 2022. [16] M.A. Sultanov, E.N. Akimova, V.E. Misilov, and Y. Nurlanuly. Parallel direct and iterative methods for solving the time-fractional diffusion equation on multicore pro- cessors. Mathematics, 10(3):323, 2022. [17] S. Kumari and M. Mehra. Numerical solution to loaded difference scheme for time- fractional diffusion equation with temporal loads. J. Math. Chem., 63(1):105–131, 2025. [18] J. Singh, A.M. Alshehri, S. Momani, S. Hadid, and D. Kumar. Computational analysis of fractional diffusion equations occurring in oil pollution. Mathematics, 10(20):3827, 2022. [19] H. Jafari, R.M. Ganji, S.M. Narsale, M. Kgarose, and V.T. Nguyen. Application of Hosoya polynomial to solve a class of time-fractional diffusion equations. Fractals, 31(04):2340059, 2023. [20] P. Roul, V.M.K.P. Goura, and R. Cavoretto. A numerical technique based on B- spline for a class of time-fractional diffusion equation. Numer. Methods Partial Differ. Equations, 39(1):45–64, 2023. [21] S. Irandoust-Pakchin, M. Hossein Derakhshan, S. Rezapour, and M. Adel. An efficient numerical method for the distributed-order time-fractional diffusion equation with the error analysis and stability properties. Math. Methods Appl. Sci., 48(3):2743–2765, 2025. [22] J.-L. Zhang, Z.-W. Fang, and H.-W. Sun. Exponential-sum-approximation tech- nique for variable-order time-fractional diffusion equations. J. Appl. Math. Comput., 48(3):323–347, 2022. [23] V.F. Morales-Delgado, J.F. Gómez-Aguilar, and M.A. Taneco-Hernandez. Analytical solution of the time fractional diffusion equation and fractional convection-diffusion equation. Rev. Mex. Fis., 65(1):82–88, 2019. [24] H. Zhu and C. Xu. A highly efficient numerical method for the time-fractional diffusion equation on unbounded domains. J. Sci. Comput., 99(2):47, 2024. [25] O. Imanuvilov, K. Ito, and M. Yamamoto. Inverse coefficient problems for one- dimensional time-fractional diffusion equations. Appl. Math. Lett., 160:109351, 2025. [26] X. Liu. Error analysis of a fully discrete method for time-fractional diffusion equations with a tempered fractional gaussian noise. J. Comput. Appl. Math., 449:115953, 2024. [27] N. Kedia, A.A. Alikhanov, and V.K. Singh. Robust finite difference scheme for the non-linear generalized time-fractional diffusion equation with non-smooth solution. Math. Comput. Simul., 219:337–354, 2024. [28] J.P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Publications, Mineola, NY, 2nd edition, 2001. [29] L.N. Trefethen. Spectral Methods in MATLAB. SIAM, Philadelphia, 2000. [30] H.M. Ahmed. Numerical solutions for singular Lane-Emden equations using shifted Chebyshev polynomials of the first kind. Contemp. Math., 4:132–149, 2023. [31] M. Abdelhakem, T. Alaa-Eldeen, D. Baleanu, M.G. Alshehri, and M. El-Kady. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 23 of 25 Approximating real-life BVPs via Chebyshev polynomials’ first derivative pseudo- Galerkin method. Fractal Fract., 5(4):165, 2021. [32] Y.H. Youssri, W.M. Abd-Elhameed, A.A. Elmasry, and A.G. Atta. An efficient Petrov–Galerkin scheme for the Euler–Bernoulli beam equation via second-kind Chebyshev polynomials. Fractal Fract., 9(2):78, 2025. [33] Y.H. Youssri, W.M. Abd-Elhameed, and M. Abdelhakem. A robust spectral treatment of a class of initial value problems using modified Chebyshev polynomials. Math. Methods Appl. Sci., 44(11):9224–9236, 2021. [34] J.C. Mason and D.C. Handscomb. Chebyshev Polynomials. Chapman and Hall, New York, NY, CRC, Boca Raton, 2003. [35] M. Masjed-Jamei. Some New Classes of Orthogonal Polynomials and Special Func- tions: A Symmetric Generalization of Sturm-Liouville Problems and its Conse- quences. PhD thesis, University of Kassel, Kassel, Germany, 2006. [36] M.M. Khader and M.M. Babatin. Evaluating the impacts of thermal conductivity on Casson fluid flow near a slippery sheet: Numerical simulation using sixth-kind Chebyshev polynomials. J. Nonlinear Math. Phys., 30(4):1834–1853, 2023. [37] K. Sadri and H. Aminikhah. A new efficient algorithm based on fifth-kind Cheby- shev polynomials for solving multi-term variable-order time-fractional diffusion-wave equation. Int. J. Comput. Math., 99(5):966–992, 2022. [38] W.M. Abd-Elhameed, Y.H. Youssri, and A.G. Atta. Adopted spectral Tau approach for the time-fractional diffusion equation via seventh-kind Chebyshev polynomials. Bound. Value Probl., 2024(1):102, 2024. [39] W.M. Abd-Elhameed, Y.H. Youssri, A.K. Amin, and A.G. Atta. Eighth-kind Cheby- shev polynomials collocation algorithm for the nonlinear time-fractional generalized Kawahara equation. Fractal Fract., 7(9):652, 2023. [40] W.M. Abd-Elhameed, O.M. Alqubori, and A.G. Atta. A collocation procedure for the numerical treatment of Fitzhugh–Nagumo equation using a kind of Chebyshev polynomials. AIMS Math., 10(1):1201–1223, 2025. [41] J. Shen, T. Tang, and L.-L. Wang. Spectral Methods: Algorithms, Analysis and Applications, volume 41 of Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg, 2011. [42] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods: Evolu- tion to Complex Geometries and Applications to Fluid Dynamics. Scientific Compu- tation. Springer, Berlin, Heidelberg, 2007. [43] J.S. Hesthaven, S. Gottlieb, and D. Gottlieb. Spectral Methods for Time-Dependent Problems, volume 21 of Cambridge Monographs on Applied and Computational Math- ematics. Cambridge University Press, Cambridge, 2007. [44] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods in Fluid Dynamics. Springer Series in Computational Physics. Springer-Verlag, Berlin, Hei- delberg, 1988. [45] S. Kumar, V. Gupta, and D. Zeidan. An efficient collocation technique based on operational matrix of fractional-order Lagrange polynomials for solving the space- time fractional-order partial differential equations. Appl. Numer. Math., 204:249–264, W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 24 of 25 2024. [46] G. Manohara and S. Kumbinarasaiah. An innovative Fibonacci wavelet collocation method for the numerical approximation of Emden-Fowler equations. Appl. Numer. Math., 201:347–369, 2024. [47] M.A. Taema, M.A. Dagher, and Y.H. Youssri. Spectral collocation method via Fer- mat polynomials for Fredholm–Volterra integral equations with singular kernels and fractional differential equations. Palestine J. Math., 14(2):481–492, 2025. [48] M.A. Taema and Y.H. Youssri. Third-kind Chebyshev spectral collocation method for solving models of two interacting biological species. Contemp. Math., 5:6189–6207, 2024. [49] A.H. Bhrawy and W.M. Abd-Elhameed. New algorithm for the numerical solutions of nonlinear third-order differential equations using Jacobi–Gauss collocation method. Math. Probl. Eng., 2011:837218, 2011. [50] W.M. Abd-Elhameed, A.M. Al-Sady, O.M. Alqubori, and A.G. Atta. Numerical treat- ment of the fractional Rayleigh-Stokes problem using some orthogonal combinations of Chebyshev polynomials. AIMS Math., 9(09):25457–25481, 2024. [51] M.A. Abbas, Y.H. Youssri, M. El-Kady, and M. Abdelhakem. High-accuracy modified spectral techniques for two-dimensional integral equations. J. Comput. Appl. Mech., 56(2):364–379, 2025. [52] S.M. Sayed, A.S. Mohamed, E.M. Abo-Eldahab, and Y.H. Youssri. Spectral frame- work using modified shifted Chebyshev polynomials of the third kind for numerical solutions of one- and two-dimensional hyperbolic telegraph equations. Bound. Value Probl., 2025(1):7, 2025. [53] R.M. Hafez, H.M. Ahmed, O.M. Alqubori, A.K. Amin, and W.M. Abd-Elhameed. Efficient spectral Galerkin and collocation approaches using telephone polynomials for solving some models of differential equations with convergence analysis. Mathematics, 13(6):918, 2025. [54] Y. H. Youssri, R. M. Hafez, and A. G. Atta. An innovative pseudo-spectral Galerkin algorithm for the time-fractional tricomi-type equation. Phys. Scr., 99(10):105238, 2024. [55] A.A. El-Sayed, S. Boulaaras, and N.H. Sweilam. Numerical solution of the fractional- order logistic equation via the first-kind Dickson polynomials and spectral Tau method. Math. Methods Appl. Sci., 46(7):8004–8017, 2023. [56] W.M. Abd-Elhameed, J.A.T. Machado, and Y.H. Youssri. Hypergeometric fractional derivatives formula of shifted Chebyshev polynomials: Tau algorithm for a type of fractional delay differential equations. Int. J. Nonlinear Sci. Numer. Simul., 23(7- 8):1253–1268, 2022. [57] W.M. Abd-Elhameed, O.M. Alqubori, A.K. Al-Harbi, M.H. Alharbi, and A.G. Atta. Generalized third-kind Chebyshev Tau approach for treating the time fractional Cable problem. Electron. Res. Archive, 32(11), 2024. [58] W.M. Abd-Elhameed, M.A. Abdelkawy, O.M. Alqubori, and A.G. Atta. An accurate tau-based spectral algorithm for the time fractional bioheat transfer model. Bound. Value Probl., 2025:124, 2025. W. M. Abd-Elhameed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6838 25 of 25 [59] Y. Xu. An integral formula for generalized Gegenbauer polynomials and Jacobi poly- nomials. Adv. Appl. Math., 29(2):328–343, 2002. [60] A. Draux, M. Sadik, and B. Moalla. Markov–Bernstein inequalities for generalized Gegenbauer weight. Appl. Numer. Math., 61(12):1301–1321, 2011. [61] Y.H. Youssri and A.G. Atta. Petrov-Galerkin Lucas polynomials procedure for the time-fractional diffusion equation. Contemp. Math., 4:230–248, 2023. [62] X. Zhao, L.L. Wang, and Z. Xie. Sharp error bounds for Jacobi expansions and Gegenbauer–Gauss quadrature of analytic functions. SIAM J. Numer. Anal., 51(3):1443–1469, 2013. [63] K. Sadri and H. Aminikhah. Chebyshev polynomials of sixth kind for solving nonlinear fractional PDEs with proportional delay and its convergence analysis. J. Funct. Spaces, 2022(1):9512048, 2022. [64] K. Sayevand, A. Yazdani, and F. Arjang. Cubic B-spline collocation method and its application for anomalous fractional diffusion equations in transport dynamic systems. J. Vib. Control, 22(9):2173–2186, 2016.