EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 4, Article Number 7095
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

A Laplace-Chebyshev Spectral Method for
Multi-Dimensional Anomalous Transport

Kamran1,∗, Bibi Zahra1, Zeeshan Ali2, Ahmad Aloqaily3, Nabil Mlaiki3
1 Department of Mathematics, Islamia College Peshawar, Peshawar 25120, Khyber Pakhtunkhwa,
Pakistan
2 Department of Information Management, National Yunlin University of Science and
Technology, Douliu, Taiwan, Republic of China
3 Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833,
Riyadh 11586, Saudi Arabia

Abstract. Anomalous transport processes, such as subsurface contaminant spread or wave at-
tenuation in viscoelastic materials, are governed by time-fractional diffusion-wave equations. The
non-local nature of fractional operators and the high computational cost of addressing multi-
dimensional spaces pose significant challenges for numerical simulations. To overcome this, we
develop a novel hybrid spectral method combining the Laplace transform (LT) technique with
the Chebyshev spectral collocation method (CSCM) for solving TFDWEs featuring the modified
Atangana-Baleanu-Caputo derivative, chosen for its non-singular kernel and efficiency in modeling
complex memory effects. Our numerical scheme, temporal and spatial discretizations, are decou-
pled. The LT handles the fractional time derivative exactly in the Laplace domain, removing
time-stepping restrictions and convolution costs, while the CSCM ensures the exponential con-
vergence in the spatial domain. The numerical inversion of LT is obtained using the improved
Talbot method, guaranteeing rapid O(e−cN) convergence. This work provides not only a robust
computational technique but also a rigorous mathematical analysis, establishing clear conditions
for the existence, uniqueness, and Ulam-Hyers stability of the solutions. The dimensional flexibility
of our technique is demonstrated through 1D, 2D, and 3D numerical examples, which confirm its
computational efficiency and high accuracy. This work provides a robust and stable numerical
approach that can be extended to model complex multi-scale transport problems across applied
mathematics and engineering.
2020 Mathematics Subject Classifications: 44A10, 65R10, 35A22
Key Words and Phrases: Diffusion-wave equation, modified Atangana-Baleanu derivative,
Laplace transform, Chebyshev spectral method, Talbot’s method, uniqueness and existence

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.7095

Email addresses: kamran.maths@icp.edu.pk (Kamran),
zeeshan@yuntech.edu.tw (Z. Ali),
maloqaily@psu.edu.sa (A. Aloqaily),
nmlaiki@psu.edu.sa; nmlaiki2012@gmail.com (N. Mlaiki)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 2 of 30

1. Introduction

Anomalous transport processes, such as subsurface pollutant dispersion or wave atten-
uation in viscoelastic materials, are common throughout scientific disciplines but cannot
be described using classical integer-order models. Time-fractional diffusion-wave equations
(TFDWEs) have emerged as an effective structure for describing these complicated phe-
nomena, effectively integrating memory effects and power-law dynamics using fractional-
order operators [1–4]. The growth of these operators is vital: from the initial work of
Abel, Riemann, and Liouville [5] to the broadly recognized Caputo derivative [6], which
added historicity but had a singular kernel. This limitation motivated the development of
non-singular alternatives, resulting in the modified Atangana-Baleanu-Caputo (MABC)
derivative [7]. The MABC derivative, with its non-singular Mittag-Leffler kernel, provides
a more reliable model for systems with complex memory and hereditary properties. The
Modified Atangana-Baleanu (MABC) derivative extends the original ABC operator to a
wider function space. Crucially, as demonstrated in [7], the MABC derivative can solve a
class of fractional differential equations intractable under the standard ABC definition. Its
applications span numerous scientific fields, including viscoelasticity [8, 9], control theory
[? ], and biological systems [10–13], with further examples available in the cited literature.

The numerical solution of TFDWEs is a dynamic and challenging area of research. A
wide range of approaches has been developed to address these problems, each with its own
advantages. Analytical solutions, such as those by Mainardi [14] for 1D cases or Agrawal
[15] for bounded domains, provide foundational insights but are impractical for complex,
multi-dimensional problems. Therefore, numerical methods have emerged, including finite
difference methods (FDM) with compact schemes [16], solvers for multi-term equations
[17], and the alternating direction implicit (ADI) method for 2D problems [18]. More
advanced techniques, such as meshless techniques [19], spectral methods [20, 21], wavelet-
based approaches [22], and hybrid Laplace-spectral methods [23], have further improved
spatial accuracy and computational efficiency. These advancements are further exemplified
by recent developments such as the Galerkin spectral method with high-order differences
[24] and fractional multi-step methods [25], underscoring the rapid progress in the field.

Despite these advances, the non-local nature of fractional operators like the MABC
poses a significant computational challenge. Time-stepping methods, including the FDM
and spectral methods, suffer significant costs in long-time or high-dimensional simulations
due to the need to store and process the entire solution history at each step. This history
dependence increases computational cost and complexity, making accurate 3D simulations
of anomalous transport processes a significant challenge. While hybrid methods, such
as those coupling the LT with a spectral collocation method [23], minimize temporal
complexity, a complete framework that fully integrates the time-fractional derivative with
rigorous mathematical analysis and demonstrates high performance in a multi-dimensional
setting has not yet been developed.

This work addresses this gap by proposing a novel hybrid algorithm that combines
the LT with the CSCM. Our approach decouples the temporal and spatial challenges
of TFDWEs: the LT transforms the MABC time-fractional operator exactly into the



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 3 of 30

Laplace domain, eliminating the convolution burden and stability constraints of time-
stepping methods. The resulting parameterized Helmholtz-type problems are solved using
the CSCM, which achieves exponential convergence in space for smooth solutions [26–28].
The time-domain solution is efficiently recovered using the improved Talbot method for
numerical inversion [29]. This approach not only reduces the computational cost but also
offers a unified framework for 1D, 2D, and 3D problems.

Beyond computational efficiency, this work provides rigorous mathematical and phys-
ical insights. First, we establish a solid theoretical foundation for the MABC-based
TFDWE by proving the existence, uniqueness, and Ulam-Hyers stability of its solutions,
thereby ensuring the model’s well-posedness a critical aspect often overlooked in numeri-
cal studies. Furthermore, by eliminating time-stepping, our method facilitates an efficient
analysis of how the fractional order α influences solution behavior in multi-dimensional
settings, offering new insights into memory effects within complex systems.

The key advantages of the proposed LT-based CSCM approach are threefold. First,
it avoids the time-step restrictions and temporal error accumulation inherent in step-wise
methods. Second, it efficiently handles the memory effects of the MABC derivative via
the Laplace transform while leveraging the CSCM for exponential convergence of smooth
spatial solutions. Finally, the method requires fewer discretization nodes for high accuracy
compared to other methods, which significantly reduces computational and storage costs
and provides a unified framework for solving 1D, 2D, and 3D problems.

The rest of this paper is organized as follows. Section 2 introduces the necessary
preliminary definitions. The analysis of existence, uniqueness, and stability is presented
in Section 3. Section 4 provides a detailed description of the proposed numerical method.
Numerical results for 1D, 2D, and 3D examples are discussed in Section 5, and finally,
conclusions are drawn in Section 6.

2. Preliminaries

This section presents the fundamental definitions and mathematical prerequisites for
the current study.

Definition 1. The Laplace transform (LT) of a function u(x̄, t) is defined as:

L {u(x̄, τ)} = û(x̄, s) =

∫ ∞

0
esτu(x̄, τ)dτ

Definition 2. The MABC derivative of order 0 < α < 1 of a function u(x̄, τ) ∈ L1(0, T )
in Caputo sense is defined by [30]:

MABC
0 Dα

τ u(x̄, τ) =
β(α)

1− α

[
u(x̄, τ)− Eα (−ηατ

α) u(x̄, 0)

− ηα

∫ τ

0
(τ − ϑ)α−1Eα,α (−ηα(τ − ϑ)α) u(x̄, ϑ)dϑ

]
,



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 4 of 30

Definition 3. Let u(n−1)(x̄, τ) ∈ L1(0, T ), then the MABC derivative of order n −
1 < α1 < n of a function in Caputo sense is defined by [30]:

MABC
0 Dα1

τ u(x̄, τ) =
β(α)

1− α

[
u(n−1)(x̄, τ)− Eα (−ηατ

α) u(n−1)(x̄, 0)

− ηα

∫ τ

0
(τ − ϑ)α−1Eα,α (−ηα(τ − ϑ)α) u(n−1)(x̄, ϑ)dϑ

]
,

where α1 = α + n − 1, ηα = α
1−α , β(α) is a normalized function having the property

β(0) = β(1) = 1, and Eα(·) is a Mittag-Leffler function defined as

Eα(τ) =
∞∑
ℓ=0

τ ℓ

Γ(αℓ+ 1)
.

Definition 4. For u(x̄, τ) ∈ L1(0, T ), the MABC fractional integral operator is defined by
[30]:

MABCT
α
0u(x̄, τ) =

1− α

β(α)
u(x̄, τ) +

α

β(α)Γ(α)

∫ τ

0
(τ − ϑ)αu(x̄, τ)dτ, τ ≥ 0.

Definition 5. The LT of MABC derivative of a function u(x̄, τ) is defined by [30]:

L
{
MABC
0 Dα

τ u(x̄, τ)
}
=

sαû(x̄, s)− sα−1û(x̄, 0)

sα(1− α) + α
.

3. Existence and uniqueness results

We establish the existence and uniqueness of the solution to the given problem in this
section. We begin by defining an appropriate functional framework: Let B(Ω,R) represent
the Banach space comprising all continuous real-valued functions defined on the compact
domain Ω = Θ× [0, 1] where Θ ⊂ R3. The compactness of Ω ensures that B(Ω,R) is well
suited for studying the continuity and boundedness properties of sections. The norm on
B(Ω,R) is a supremum norm given by:

∥u∥∞ = sup{|u(x̄, τ)| : (x̄, τ) ∈ Ω},

where x̄ ∈ Θ represents the spatial co ordinates and τ ∈ [0, 1] denotes the temporal
variable. We study the fractional diffusion-wave equation of the form [31]:

MABC
0 Dα

τ u(x̄, τ) = λ1∇2u(x̄, τ)− λ2u(x̄, τ) + f(x̄, τ), 1 < α ≤ 2, x̄ ∈ Θ, (1)

with boundary conditions
Bu(x̄, τ) = g1(x̄, τ), x̄ ∈ ∂Θ, (2)

and initial conditions
u(x̄, 0) = g2(x̄),

uτ (x̄, 0) = g3(x̄),
(3)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 5 of 30

where ∇2 is the Laplacian operator, defined as: for 1D problem, x̄ = x and ∇2 = ∂2

∂x2 ; for
a 2D problem, x̄ = (x, y) and ∇2 = ∂2

∂x2 + ∂2

∂y2
; and for a 3D problem, x̄ = (x, y, z) and

∇2 = ∂2

∂x2 +
∂2

∂y2
+ ∂2

∂z2
. Here Θ is the domain and ∂Θ is its boundary. The constants λ1 and

λ2 are arbitrary, and the forcing term f(x̄, τ) is sufficiently smooth. The functions g1(x̄, τ),
g2(x̄), and g3(x̄) are given continuous functions. B is the boundary differential operator,
and MABC

0 Dα
τ u(x̄, τ), denotes the MABC derivative of order 1 < α < 2. Applying the

MABC fractional integral operator to Eq (1), yields:

u(x̄, τ) = g2(x̄) + g3(x̄)t+
(1− α)

β(α)

(
λ1∇2u(x̄, τ)− λ2u(x̄, τ) + f(x̄, τ)

)
+

α

β(α)Γ(α)

(∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ)− λ2u(x̄, ϑ) + f(x̄, ϑ)

)
dϑ

)
.

(4)

Next, we define the operator J : B(Ω,R) → B(Ω,R), transforming the problem into a
fixed-point formulation:

J u(x̄, τ) = g2(x̄) + g3(x̄)t+
(1− α)

β(α)

(
λ1∇2u(x̄, τ)− λ2u(x̄, τ) + f(x̄, τ)

)
+

α

β(α)Γ(α)

(∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ)− λ2u(x̄, ϑ) + f(x̄, ϑ)

)
dϑ

)
.

(5)

The fixed point of the operator corresponds to the solution to Eqs. (1)–(3). We introduce
the following hypotheses for any (ϑ, τ) ∈ (Ω):

(H1) |∇2u(x̄, τ)| ≤ ϵ1|u(x̄, τ)|,
(H2) |g2| ≤ ϵ2,

(H3) |g3| ≤ ϵ3,

(H4) |f(x̄, τ)| ≤ ϵ4,

(H5) |∇2u1(x̄, τ)−∇2u2(x̄, τ)| ≤ ϵ5|u1(x̄, τ)− u2(x̄, τ)|,
(H6) |∇2u1(x̄1, τ1)−∇2u2(x̄2, τ2)| ≤ Lf1 (|x̄1 − x̄2|+ |τ1 − τ2||) ,
(H7) |u1(x̄1, τ1)− u2(x̄2, τ2)| ≤ Lf2 (|x̄1 − x̄2|+ |τ1 − τ2||) ,
(H8) |f(x̄1, τ1)− f(x̄2, τ2)| ≤ Lf3 (|x̄1 − x̄2|+ |τ1 − τ2||) ,
(H9) |g2(x̄1)− g2(x̄2)| ≤ Lg2 |x̄1 − x̄2| ≤ Lg2δ1,

(H10) |g3(x̄1)− g3(x̄2)| ≤ Lg3 |x̄1 − x̄2| ≤ Lg3δ1,

where ϵ1, ϵ2, ϵ3, ϵ4, ϵ5,Lf1 ,Lf2 ,Lf3 ,Lg2 ,Lg3 > 0 are constants.
The problem defined by Eqs. (1)–(3) has at least one solution.
Proof. The proof proceeds in several steps. We apply the Schaefer fixed-point theorem

to establish the existence of a solution.
Step 1: In the first we show that the operator J is continuous. Consider a sequence



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 6 of 30

um → u in B(Ω,R). For (x̄, τ) ∈ Ω, we compute:

∥J um(x̄, τ)− J u(x̄, τ)∥∞ = sup
(x̄,τ)∈Ω

{
|J um(x̄, τ)− J u(x̄, τ)|

}
= sup

(x̄,τ)∈Ω

{∣∣∣∣1− α

β(α)

(
λ1∇2um(x̄, τ)− λ2um(x̄, τ)

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
λ1∇2um(x̄, ϑ)− λ2um(x̄, ϑ)

)
dϑ

− 1− α

β(α)

(
λ1∇2u(x̄, τ)− λ2u(x̄, τ)

)
− α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ)− λ2u(x̄, ϑ)

)
dϑ

∣∣∣∣}
≤ sup

(x̄,τ)∈Ω

{
1− α

β(α)

(
|λ1||∇2um(x̄, τ)−∇2u(x̄, τ)|+ |λ2||um(x̄, τ)− u(x̄, τ)|

)
− α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
|λ1||∇2um(x̄, ϑ)−∇2u(x̄, ϑ)|+ |λ2||um(x̄, ϑ)− u(x̄, ϑ)|

)
dϑ.

}
Using hypotheses H5, we have

∥J um(x̄, τ)− J u(x̄, τ)∥∞ ≤ sup
(x̄,τ)∈Ω

{
1− α

β(α)

(
|λ1|ϵ5|um(x̄, τ)− u(x̄, τ)|+ |λ2||um(x̄, τ)− u(x̄, τ)|

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
ϵ5|λ1||um(x̄, ϑ)− u(x̄, ϑ)|+ |λ2||um(x̄, ϑ)− u(x̄, ϑ)|

)
dϑ

}
≤1− α

β(α)

(
(|λ1|ϵ5 + |λ2|)∥um − u∥∞

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
(ϵ5|λ1|+ |λ2|)∥um − u∥∞

)
dϑ

=
1− α

β(α)

(
(|λ1|ϵ5 + |λ2|)∥um − u∥∞

)
+

τα

Γ(α)β(α)

(
(ϵ5|λ1|+ |λ2|)∥um − u∥∞

)
=
(1− α)Γ(α) + τα

Γ(α)β(α)

(
(|λ1|ϵ5 + |λ2|)∥um − u∥∞

)
Since um → u in B(Ω,R), we have∥um − u∥∞ → 0, as m → ∞. Hence, ∥J um −
J u∥∞ → 0, proving that J is continuous.
Step 2: Boundedness of J . Let Rγ = {u ∈ C(Λ,R) : ∥u∥∞ ≤ γ} for some γ > 0. For
u ∈ Rγ , we estimate:

|J u(x̄, τ)| =
∣∣∣∣g2(x̄) + g3(x̄)τ +

1− α

β(α)

(
λ1∇2u(x̄, τ)− λ2u(x̄, τ) + f(x̄, τ)

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ)− λ2u(x̄, ϑ) + f(x̄, ϑ)

)
dϑ

∣∣∣∣
≤|g2(x̄)|+ |g3(x̄)|τ +

1− α

β(α)

(
|λ1||∇2u(x̄, τ)|+ |λ2||u(x̄, τ)|+ |f(x̄, τ)|

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
|λ1||∇2u(x̄, ϑ)|+ |λ2||u(x̄, ϑ)|+ |f(x̄, ϑ)|

)
dϑ

∣∣∣∣.
Using H1 −H4, we have



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 7 of 30

|J u(x̄, τ)| ≤ ϵ2 + τϵ3 +
1− α

β(α)

(
|λ1|ϵ1|u(x̄, τ)|+ |λ2||u(x̄, τ)|+ ϵ4

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
|λ1|ϵ1|u(x̄, ϑ)|+ |λ2||u(x̄, ϑ)|+ ϵ4

)
dϑ

∣∣∣∣
or

∥J u∥∞ ≤ϵ2 + τϵ3 +
1− α

β(α)

(
|λ1|ϵ1∥u∥∞ + |λ2|∥u∥∞ + ϵ4

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
|λ1|ϵ1∥u∥∞ + |λ2|∥u∥∞ + ϵ4

)
dϑ

∣∣∣∣
= ϵ2 + τϵ3 +

Γ(α)(1− α) + τα

Γ(α)β(α)

(
|λ1|ϵ1∥u∥∞ + |λ2|∥u∥∞ + ϵ4

)
≤ ϵ2 + τϵ3 +

Γ(α)(1− α) + τα

Γ(α)β(α)

(
|λ1|ϵ1γ + |λ2|γ + ϵ4

)
,

since τ ∈ [0, 1], we have

∥J ≤ u∥∞ ≤ ϵ2 + ϵ3 +
Γ(α)(1− α) + 1

Γ(α)β(α)

(
|λ1|ϵ1γ + |λ2|γ + ϵ4

)
=: ρab,

where ρab is a constant independent of u. Thus, J is bounded.

Step 3: Equicontinuity of J . Let u ∈ Rγ and (x̄1, τ1), (x̄2, τ2) ∈ Ω with ∥x̄1−x̄2∥ < δ1
and |τ1 − τ2| < δ2. We analyze:

|J u(x̄1, τ1)− J u(x̄2, τ2)| ≤ |g2(x̄1)− g2(x̄2)|+ |g3(x̄1)− g3(x̄2)|τ1 + |g3(x̄2)||τ1 − τ2|

+
1− α

β(α)

∣∣∣∣λ1∇2u(x̄1, τ1)− λ1∇2u(x̄2, τ2) + λ2u(x̄2, τ2)− λ2u(x̄1, τ1) + f(x̄1, τ1)− f(x̄2, τ2)

∣∣∣∣
+

α

Γ(α)β(α)

∣∣∣∣ ∫ τ1

0
(τ1 − ϑ)α−1

(
λ1∇2u(x̄1, ϑ)− λ2u(x̄1, ϑ) + f(x̄1, ϑ)

)
dϑ

−
∫ τ2

0
(τ2 − ϑ)α−1

(
λ1∇2u(x̄2, ϑ)− λ2u(x̄2, ϑ) + f(x̄2, ϑ)

)
dϑ

∣∣∣∣
Using H3 and H6 −H10, we have

|J u(x̄1, τ1)− J u(x̄2, τ2)| ≤ Lg2δ1 + Lg3δ1 + ϵ3δ2 +
1− α

β(α)
(|λ1|Lf1 (|x̄1 − x̄2|+ |τ1 − τ2||) + |λ2|Lf2 (|x̄1 − x̄2|

+|τ1 − τ2||) + Lf3 (|x̄1 − x̄2|+ |τ1 − τ2||))

+
α

Γ(α)β(α)

∣∣∣∣ ∫ τ1

0
(τ1 − ϑ)α−1

(
λ1∇2u(x̄1, ϑ)− λ2u(x̄1, ϑ) + f(x̄1, ϑ)

)
dϑ

−
∫ τ2

0
(τ2 − ϑ)α−1

(
λ1∇2u(x̄2, ϑ)− λ2u(x̄2, ϑ) + f(x̄2, ϑ)

)
dϑ

∣∣∣∣,



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 8 of 30

Next, we split the difference of integrals into three parts as follows:∣∣∣∣ ∫ τ1

0
(τ1 − ϑ)α−1

(
λ1∇2u(x̄1, ϑ)− λ2u(x̄1, ϑ) + f(x̄1, ϑ)

)
dϑ

−
∫ τ2

0
(τ2 − ϑ)α−1

(
λ1∇2u(x̄2, ϑ)− λ2u(x̄2, ϑ) + f(x̄2, ϑ)

)
dϑ

∣∣∣∣
=

[ ∫ τ1

0

[
(τ1 − ϑ)α−1 − (τ2 − ϑ)α−1

] (
λ1∇2u(x̄1, ϑ)− λ2u(x̄1, ϑ) + f(x̄1, ϑ)

)
dϑ︸ ︷︷ ︸

(i)

+

∫ τ1

0
(τ2 − ϑ)α−1

(
λ1[∇2u(x̄1, ϑ)−∇2u(x̄2, ϑ)]− λ2[u(x̄1, ϑ)− u(x̄2, ϑ)] + [f(x̄1, ϑ)− f(x̄2, ϑ)]

)
dϑ︸ ︷︷ ︸

(ii)

+

∫ τ2

τ1

(τ2 − ϑ)α−1
(
λ1∇2u(x̄2, ϑ)− λ2u(x̄2, ϑ) + f(x̄2, ϑ)

)
dϑ︸ ︷︷ ︸

(iii)

]

for the integral (i), using the mean value theorem, we have

|(τ1 − ϑ)α−1 − (τ2 − ϑ)α−1| ≤ (α− 1)(τ2 − τ1)(τ1 − ϑ)α−2

Therefore:

(i) ≤ (λ1∥∇2u∥∞ + λ2∥u∥∞ + ∥f∥∞)(α− 1)δ2

∫ τ1

0
(τ1 − ϑ)α−2dϑ

= (λ1∥∇2u∥∞ + λ2∥u∥∞ + ∥f∥∞)δ2τ
α−1
1 .

Similarly, using Lipschitz conditions for integral (ii), we have:

(ii) ≤ (λ1Lf1 + λ2Lf2 + Lf3)δ1

∫ τ1

0
(τ2 − ϑ)α−1dϑ

= (λ1Lf1 + λ2Lf2 + Lf3)δ1
(τα2
α

− δ2
α

)
.

and using direct estimation, we have:

(iii) ≤ (λ1∥∇2u∥∞ + λ2∥u∥∞ + ∥f∥∞)

∫ τ2

τ1

(τ2 − ϑ)α−1dϑ

= (λ1∥∇2u∥∞ + λ2∥u∥∞ + ∥f∥∞)
δ2
α
.

Combining all the results, we obtain:



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 9 of 30

∥J u(x̄1, τ1)− J u(x̄2, τ2)∥∞ ≤ Lg2δ1 + Lg3δ1 + ϵ3δ2 +
1− α

β(α)

(
|λ1|Lf1(δ1 + δ2) + |λ2|Lf2(δ1 + δ2)

+ Lf3(δ1 + δ2)

)
+

α

Γ(α)β(α)

(
(λ1∥∇2u∥∞ + λ2∥u∥∞ + ∥f∥∞)δ2τ

α−1
1 + (λ1Lf1 + λ2Lf2 + Lf3)δ1

(τα2
α

− δ2
α

)
+ (λ1∥∇2u∥∞ + λ2∥u∥∞ + ∥f∥∞)

δ2
α

)
.

Hence
∥J u(x̄1, τ1)− J u(x̄2, τ2)∥∞ → 0

as δ1, δ2 → 0, proving equicontinuity.
Step 4: A Priori Bound. Define ℵ = {u ∈ B(Ω,R) : u = εJ u, ε ∈ (0, 1)}. For u ∈ ℵ,

we have
|u| = |εJ u| = ε|J u| ≤ ερab,

where ρab is defined as in Step 2. The inequality ∥u∥∞ ≤ ρab implies the operator ℵ is
bounded. Therefore, by the Schaefer fixed-point theorem [32], J has at least one fixed
point, ensuring the existence of a solution to the problem.

The problem defined in Eq. (1) has a unique solution if the following condition is
satisfied:

(1− α)Γ(α) + τα

Γ(α)β(α)

((
ϵ5|λ1|+ |λ2|

)
∥u1 − u2∥∞

)
< 1. (6)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 10 of 30

Proof.

∥J u1(x̄, τ)− J u2(x̄, τ)∥∞ = sup

{
|J u1(x̄, τ)− J u2(x̄, τ)|

}
= sup

{∣∣∣∣1− σ

β(α)

(
λ1∇2u1(x̄, τ)− λ2u1(x̄, τ)

)
+

σ

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u1(x̄, ϑ)− λ2u1(x̄, ϑ)

)
dϑ

− 1− α

β(α)

(
λ1∇2u2(x̄, ϑ)− λ2u2(x̄, ϑ)

)
− α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u2(x̄, ϑ)− λ2u2(x̄, ϑ)

)
dϑ

∣∣∣∣}
≤ sup

{
1− α

β(α)

(
|λ1||∇2u1(x̄, τ)−∇2u2(x̄, τ)|+ |λ2||u1(x̄, τ)− u2(x̄, τ)|

)
− α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
|λ1||∇2u1(x̄, ϑ)−∇2u2(x̄, ϑ)|+ |λ2||u1(x̄, ϑ)− u2(x̄, ϑ)|

)
dϑ

}
≤ sup

{
1− α

β(α)

(
ϵ5|λ1||u1(x̄, τ)− u2(x̄, τ)|+ |λ2||u1(x̄, τ)− u2(x̄, τ)|

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

(
ϵ5|λ1||u1(x̄, τ)− u2(x̄, τ)|+ |λ2||u1(x̄, τ)− u2(x̄, τ)|

)
dϑ

}
=

1− α

β(α)

((
ϵ5|λ1|+ |λ2|

)
∥u1 − u2∥∞

)
+

α

Γ(α)β(α)

∫ τ

0
(τ − ϑ)α−1

((
ϵ5|λ1|+ |λ2|

)
∥u1 − u2∥∞

)
dϑ

=
(1− α)Γ(α) + τα

Γ(α)β(α)

((
ϵ5|λ1|+ |λ2|

)
∥u1 − u2∥∞

)
.

If condition (6) holds, then the operator J is a contraction. By the Banach fixed-point
theorem [32], this implies that J admits a unique fixed point, which in turn implies that
problem (1) has a unique solution.

4. Ulam-Hyers stability

The solution of problem (1)–(3) is Ulam-Hyers stable under hypotheses H5.
Proof. Let the exact solution of the problem is given by:

u(x̄, τ) = g2(x̄) + g3(x̄)t+
(1− α)

β(α)

(
λ1∇2u(x̄, τ) + λ2u(x̄, τ) + f(x̄, τ)

)
+

α

β(α)Γ(α)

(∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ) + λ2u(x̄, ϑ) + f(x̄, ϑ)

)
dϑ

)
,

(7)

and let u(x̄, τ) be the approximate solution defined as:

u(x̄, τ) = g2(x̄) + g3(x̄)t+
(1− α)

β(α)

(
λ1∇2u(x̄, τ) + λ2u(x̄, τ) + f(x̄, τ)

)
+

α

β(α)Γ(α)

(∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ) + λ2u(x̄, ϑ) + f(x̄, ϑ) + F(x̄, ϑ)

)
dϑ

)
,

(8)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 11 of 30

where F(x̄, ϑ) is a perturbation term with |F(x̄, ϑ)| ≤ ϵ7 for some constant ϵ7 > 0.
Subtracting (8) from (7), we get:

|u(x̄, τ)− u(x̄, τ)| =
∣∣∣∣(g2(x̄) + g3(x̄)t+

(1− α)

β(α)

(
λ1∇2u(x̄, τ) + λ2u(x̄, τ) + f(x̄, τ)

)
+

α

β(α)Γ(α)

(∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ) + λ2u(x̄, ϑ) + f(x̄, ϑ)

)
dϑ

))
−
(
g2(x̄) + g3(x̄)t+

(1− α)

β(α)

(
λ1∇2u(x̄, τ) + λ2u(x̄, τ) + f(x̄, τ)

)
+

α

β(α)Γ(α)

(∫ τ

0
(τ − ϑ)α−1

(
λ1∇2u(x̄, ϑ) + λ2u(x̄, ϑ) + f(x̄, ϑ) + F(x̄, ϑ)

)
dϑ

))∣∣∣∣
≤ (1− α)

β(α)

(
|λ1||∇2u(x̄, τ)−∇2u(x̄, τ)|+ |λ2||u(x̄, τ)− u(x̄, τ)|

)
+

α

β(α)Γ(α)

∫ τ

0
(τ − ϑ)α−1

(
|λ1||∇2u(x̄, τ)−∇2u(x̄, τ)|+ |λ2||u(x̄, τ)− u(x̄, τ)|

+|F(x̄, τ)|) dϑ

≤ (1− α)

β(α)
(|λ1|ϵ5|u(x̄, τ)− u(x̄, τ)|+ |λ2||u(x̄, τ)− u(x̄, τ)|)

+
α

β(α)Γ(α)

∫ τ

0
(τ − ϑ)α−1 (|λ1|ϵ5|u(x̄, τ)− u(x̄, τ)|+ |λ2||u(x̄, τ)− u(x̄, τ)|+ ϵ7) dϑ

taking the supremum norm:

∥u− u∥∞ ≤ (1− α)

β(α)
(|λ1|ϵ5∥u− u∥∞ + |λ2|∥u− u∥∞)

+
α

β(α)Γ(α)

∫ τ

0
(τ − ϑ)α−1 (|λ1|ϵ5∥u− u∥∞ + |λ2|∥u− u∥∞ + ϵ7) dϑ

≤
(
(1− α)

β(α)
+

τα

β(α)Γ(α)

)
(|λ1|ϵ5∥u− u∥∞ + |λ2|∥u− u∥∞)

+
ϵ7τ

α

β(α)Γ(α)

)
=

(
((1− α)Γ(α) + τα)(|λ1|ϵ5 + |λ2|)

β(α)Γ(α)

)
∥u− u∥∞ +

ϵ7τ
α

β(α)Γ(α)
.

Rearranging, we get:[
1−

(
((1− α)Γ(α) + τα)(|λ1|ϵ5 + |λ2|)

β(α)Γ(α)

)]
∥u− u∥∞ ≤ ϵ7τ

α

β(α)Γ(α)
.

Since Eq.(6) ensures that (1−α)Γ(α)+τα

Γ(α)β(α)

((
ϵ5|λ1|+ |λ2|

)
∥u1 − u2∥∞

)
< 1, we have

∥u− u∥∞ ≤ Gϵ.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 12 of 30

where
G =

τα{
β(α)Γ(α)− [((1− α)Γ(α) + τα)(|λ1|ϵ5 + |λ2|)]

} .

Since G > 0, the solution is Ulam-Hyers stable.

5. Proposed numerical method

The proposed numerical method for solving the TFDWE with the MABC derivative
consists of three main steps: (a) discretization of the time variable via the Laplace trans-
form, which converts the problem to the Laplace domain; (b) solution of the resulting
boundary value problem in the Laplace domain using the CSCM; and (c) recovery of
the time-domain solution by applying a numerical inverse Laplace transform based on a
modified Talbot contour and the midpoint rule.

5.1. Laplace transform
The LT is applied to TFDWE (1)–(3) to discretize the time variable. Applying the LT

to the model, we have:

β(α)
(
sαû(x̄, s)− sα−1u(x̄, 0)− sα−2ut(x̄, 0)

)
sα(1− α) + α

− λ1∇2û(x̄, s)− λ2û(x̄, s) = f̂(x̄, s), x̄ ∈ Θ

for x̄ ∈ Θ with boundary condtions:

Bû(x̄, s) = ĝ1(x̄, s), x̄ ∈ ∂Θ,

the above expression can be written in operator form as:{(
β(α)sα

sα(1− α) + α

)
I − λ1£− λ2I

}
û(x̄, s) = Ĥ(x̄, s), (9)

where:
Ĥ(x̄, s) =

β(α)sα−1g2(x̄)

sα(1− α) + α
+

β(α)sα−2g3(x̄)

sα(1− α) + α
+ f̂(x̄, s)

and the boundary conditions remain

Bû(x̄, s) = ĝ1(x̄, s), x̄ ∈ ∂Θ. (10)

Here, L = ∇2 denotes the Laplacian operator. The spatial operators in Eqs. (9)–
(10) are then discretized using the CSCM, which transforms the problem into a system of
linear equations in the Laplace domain. This system is solved for each value of the Laplace
parameter s. Finally, the solution in the time domain, u(x̄, t), is recovered by applying a
numerical inverse Laplace transform.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 13 of 30

5.2. Spectral method
The CSCM is employed in this section to discretize the spatial operators in the trans-

formed system given in Eqs. (9)–(10). This method uses Lagrange interpolation polyno-
mials (LIPs) based on Chebyshev nodes to approximate the solution over the domain Θ,
which is [−1, 1] in 1D, [−1, 1]2 in 2D, and [−1, 1]3 in 3D, with x̄ = x, (x, y), or (x, y, z)
respectively. For the 1D case the solution û(x, s) is approximated as [33, 34]:

IN(x) =

N∑
=0

ℓ(x)û(x, s),

where ℓ(x) are LIPs defined as Chebyshev nodes x, which are given by:

x =

{
cos

(
π

N

)}N

=0

. (11)

and

ℓ(x) =
N∏

j=0,j ̸=

x− xj
x − xj

, (12)

The derivative ∂û(x)
∂x is approximated using the differentiation matrix N, with elements:

{N}k, = ℓ′(xk), k, = 0, 1, 2, ...,N,

where the off-diagonal entries are:

{N}k, =
υ

υk(xk − x)
, k ̸= ,

with υ−1 =
∏N

k=0,k ̸=(xk − x), and the diagonal entries:

{N}k, = −
N∑

=0, ̸=k

{N}k,, k = 0, 1, 2, ..,N.

Higher derivatives are obtained as:

{(m)
N }k, = ℓ(m)(xk).

For 1D, £Disc = ∂2

∂x2 , approximated by 2
N. The authors of [35, 36] developed an effective

and precise formulation for constructing differentiation matrices. Specifically, [35] provides
a practical method for deriving the matrix m

N , presented as follows:

{(m)
N }k =

m

xk − x

(
υ

υk
{(m−1)
N }kk − {(m−1)

N }k
)
, k ̸= .



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 14 of 30

For Γ = [−1, 1]2, the points x̄k are presented as follows:

x̄k =

(
cos

(
π

N

)
, cos

(
πk

N

))
, , k = 0, 1, 2, ...,N.

The LIPs are:
ℓk(x̄) = ℓ(x)ℓk(y), , k = 0, 1, 2, ...,N. (13)

where ℓk(x̄k) = σk. The 2nd-order derivatives of the LIPs (13) are given as:

∂2ℓk(x̄rp)

∂x2
= ℓ′′(xr)ℓk(yp) = {2N}rσkp,

∂2ℓk(x̄rp)

∂y2
= ℓ(xr)ℓ

′′
k(yp) = σr{2N}pk,

where 2
N is the 2nd order differentiation matrix. Employing £ on ℓk(x̄rp) at x̄rp gives

£(ℓk(x̄rp)) =

(
{2N}rσkp + σr{2N}pk

)
(14)

Thus, the discretized representation of the linear differential operator £, derived using the
CSCM, is given as:

£Disc = IN ⊗ 2
N + 2

N ⊗ IN, (2D) (15)

where ⊗ represent the Kronecker product. In 3D, the domain is Γ = [−1, 1]3, and the
Chebyshev nodes are:

x̄ℓkm =

(
cos

(
πℓ

n

)
, cos

(
πk

n

)
, cos

(πm
n

))
, ℓ, k,m = 0, 1, . . . , n.

The associated LIPs in 3D are:

ℓkm(x̄) = ℓ(x)ℓk(y)ℓm(z), , k,m = 0, 1, . . . , n, (16)

where ℓkm(x̄km) = σkm, and the second-order derivatives are obtained as:

∂2ℓkm(x̄rpτ )

∂x2
= ℓ′′(xr)ℓk(yp)ℓm(zτ ) = {2n}rσkpσmt,

∂2ℓkm(x̄rpτ )

∂y2
= ℓ(xr)ℓ

′′
k(yp)ℓm(zτ ) = σr{2n}pkσmt,

∂2ℓkm(x̄rpτ )

∂z2
= ℓ(xr)ℓk(yp)ℓ

′′
m(zτ ) = σrσkp{2n}tm,

where x̄rpτ = (xr, yp, zτ ). Applying £ to ℓkm(x̄rpτ ) at x̄rpτ gives:

£(ℓkm(x̄rpτ )) =
(
{2n}rσkpσmτ + σr{2n}pkσmτ + σrσkp{2n}τm

)
, (17)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 15 of 30

The discrete 3D Laplacian is

£Disc = IN ⊗ IN ⊗ 2
N + IN ⊗ 2

N ⊗ IN + 2
N ⊗ IN ⊗ IN, (3D), (18)

using matrix £Disc in Eq. (9), we obtain the discretized system as:{(
β(α)sα

sα(1− α) + α

)
I − λ1£Disc − λ2I

}
û(x̄, s) = Ĥ(x̄, s), (19)

The conditions in Eq (10) are incorporating by considering the interpolation matrix £Disc

and considering all points x̄. Furthermore, the rows of £Disc in correspondence with
boundary nodes are replaced with unit vectors that have a one in accordance with the
diagonal elements of £Disc. Hence, the boundary conditions Bû(x̄, s) = ĝ1(x̄, s) in Eq
(10) will be implemented directly [33]. Rearranging the columns and rows of the matrix
£Disc, the following block matrix is obtained.

£Γ =

[
W F
0 I

]
,

where the non-zero block W and I is of size having order (N−NB)×(N−NB) and NB×NB.
Here NB denotes the boundary nodes. Thus, the system (9)-(10) has the following form:

£Γû(x̄, s) =

[
Ĥ(x̄, s)

f̂(x̄, s)

]
. (20)

The solution û(x̄, s) in the Laplace domain is determined by solving (20). The solution
u(x̄, t) of the original problem (1)–(3) is then recovered by applying the inverse Laplace
transform to û(x̄, s) as follows:

u(x̄, τ) =
1

2πi

∫ µ+i∞

µ−i∞
esτ û(x̄, s)ds =

1

2πi

∫
Γc

esτ û(x̄, s)ds, Re(s) > 0, (21)

5.3. Contour selection and quadrature
One of the most effective approaches for computing Eq. (21) is to deform the inte-

gration contour Γc into the left half plane to ensure the integrand decays, followed by
applying the quadrature methods. This idea traces back to 1950, originating in the work
of Talbot’s doctoral student Green [37]. Talbot later published a paper [38] in which he
generalized and improved the work of Green. Talbot suggested to deform the integration
contour into a contour Γc that begins and ends in the left plane, such that Re(s) −→ −∞
as |s| → ∞. In literature, many popular contours have been proposed, such as Talbot’s
contour [29, 38], the parabolic and hyperbolic contours [39]. The work utilizes the modified
Talbot’s contour proposed in [29] as:

Γc : s = s(ξ), ξ ∈ [−π, π], Re
(
s(±π)

)
= −∞, (22)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 16 of 30

we have
s(ξ) =

MT

τ
ε(ξ), ε(ξ) = −θ1 + θ2ξ cot(θ3ξ) + θ4iξ, (23)

where the user will select the parameters θ1, θ2, θ3, and θ4. From (23) and (21), we have

u(x̄, τ) =
1

2πi

∫ π

−π
es(ξ)τ û(x̄, s(ξ))s′(ξ)dξ. (24)

Fast and accurate approximation of (24) can be achieved using trapezoidal or midpoint
rule [29]. This work focuses on the midpoint rule with a step ℏ = 2π

MT
given by:

uapp(x̄, τ) ≈
1

MTi

MT∑
k=1

es(ξk)τ û(x̄, s(ξk))s
′(ξk), ξk = −π + (

2k − 1

2
)ℏ. (25)

5.3.1. Error analysis

The error analysis is performed in three stages:
Step 1: In the first step, the LT is employed, which transforms the given problem into
a time-independent problem, since this transformation is exact, no error is introduced in
this step.
Step 2: In step 2, we employ the CSCM discretization technique to solve the transformed
problem, with the corresponding error estimate developed below:
Utilizing the points in Eq. (11) and the LPs in Eq. (12), the interpolation operator
presented in [34] is expressed as:

IN : C(Θ) → PN, IN(û) =

N∑
=0

û(x, s)ℓ(x). (26)

Following the technique in [40], we establish the error bound. Let QN be a constant; then,
the stability estimate is expressed as:

∥IN(û)∥∞ ≤ QN∥û∥∞, ∀ û ∈ C[−1, 1]. (27)

Additionally,
IN(û) = û, for all û ∈ PN. (28)

For Chebyshev interpolation, the stability constant grows logarithmically with N:

QN = 1 +

(
ln(1 + N)π

2

)
≤ (N + 1). (29)

For any û ∈ CN+1[−1, 1], the interpolation error bound is expressed as [? ]:

∥û− IN(û)∥∞ ≤ 2−N

Γ(N + 2)
∥ûN+1∥∞. (30)



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 17 of 30

Theorem 1. [40] If û ∈ C(N+1)[−1, 1], then for q = 0, 1, ...,N

∥û(q) − IN(û)
(q)∥∞ ≤

2(Q
(q)
N + 1)

Γ(N− q + 2)

(
1

2

)(N−q+1)

∥û(N+1)∥∞, (31)

where
Q

(q)
N =

QN

Γ(q + 1)

(
Γ(N + 1)

Γ(N− q + 1)

)
.

Application to 1D case:

For 1D operator £û = ∂2û(x̄,t)
∂x2 , the error bound is expressed as:

E =

∥∥∥∥{( β(α)sα

sα(1− α) + α

)
I − λ1£− λ2I

}
û−

{(
β(α)sα

sα(1− α) + α

)
I − λ1£− λ2I

}
INû

∥∥∥∥
∞

=

∥∥∥∥( β(α)sα

sα(1− α) + α

)
(û− INû)− λ1£(û− INû)− λ2(û− INû)

∥∥∥∥
∞

≤
∣∣∣∣( β(α)sα

sα(1− α) + α

)∣∣∣∣∥û− INû∥∞ + |λ1|∥£(û− INû)∥∞ + |λ2|∥û− INû∥∞,

since sα, and β(α) are constants, so we have

E ≤
((

β(α)sα

sα(1− α) + α

)
+ λ2

)
2−N

Γ(N + 2)
∥ûN+1∥∞ + |λ1|

2(Q
(2)
N + 1)

Γ(N)

(
1

2

)(N−1)

∥û(N+1)∥∞.

Hence, we obtain
E ≤ K∥û(N+1)∥∞,

where K is a constant combining all coefficients of ∥û(N+1)∥∞. For higher dimensions
see [40].

Step 3: The final step is the numerical approximation of the integral in Eq.(24) via
Talbot’s method, implemented with the midpoint rule, whose convergence rate depends
on the following factors:

• The contour selection,

• The step size.

Parameters for optimal accuracy as determined in [29] are:

θ1 = 0.61220, θ2 = 0.50170, θ3 = 0.640700, and θ4 = 0.26450s,

and corresponding error estimate:

Eest = |uapp(x̄, τ)− u(x̄, τ)| = O(exp((−1.35800)MT)).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 18 of 30

6. Numerical Experiments

The performance of the proposed method is evaluated using three numerical examples.
Accuracy is measured with two error metrics: the absolute error (LAbs) and the maximum
absolute error (L∞), defined as:

LAbs =

∣∣∣∣u(x̄k, τ)− uapp(x̄k, τ)

∣∣∣∣,
L∞ = max

1≤k≤N

∣∣∣∣u(x̄k, τ)− uapp(x̄k, τ)

∣∣∣∣,
where u(x̄, τ) and uapp(x̄, τ) denote the exact and approximate solutions, respectively.
Here, N and MQ denote the number of Chebyshev nodes and quadrature nodes, respec-
tively. All simulations employ fixed parameters λ1 = 1 and λ2 = 0. For each example,
the source term, initial conditions, and boundary conditions are derived from the exact
solution.

Example 1

In the first example, we consider the 1D version of (1)–(3) with λ1 = 1, λ2 = 0, and
exact solution u(x, τ) = sin(πx)τ2. The performance of the proposed numerical method is
evaluated through comprehensive error analysis and computational benchmarks. Table 1
shows the L∞ error norms for varying Chebyshev nodes N and quadrature points (MQ),
demonstrating both computational efficiency and high accuracy. The solution accuracy
is verified in Figure 1a, where the upper panel shows excellent agreement between exact
and approximate solutions, while the lower panel depicts the corresponding pointwise
absolute error distribution, confirming the method’s high accuracy. Further, Figure 1b
further examines solution accuracy for different values of α, showing consistent and stable
numerical behavior.

Two major features are shown by convergence analysis: Up to MQ = 36, Figure 2a
demonstrates ideal quadrature convergence, after that, there is a slight increase in the
error, most likely as a result of numerical conditioning effects. Similarly, spectral con-
vergence with respect to spatial discretization is shown in Figure 2b, with slight error
variation at higher N values caused by the round-off errors in the Chebyshev differentia-
tion matrices.

Figures 3a and 3b quantify parametric sensitivity by plotting the L∞ error dependence
on t and α, respectively, both showing accurate results. The surface and contour plots
(Figures 4a and 4b) in the α−τ plane show the whole error nature and provide a thorough
proof of the method’s stability over the whole domain. Overall, these findings demonstrate
that the proposed scheme achieves the following: (i) uniform stability for various values of
α; (ii) exponential convergence in quadrature approximation; and (iii) spectral accuracy
in spatial discretization. Thus, the method represents an effective technique for Solving
fractional-order PDE problems.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 19 of 30

Table 1: Errors norms for Example 1 with varying MQ, α, and N.

MQ N
α = 1.5 α = 1.75

L∞ C.Time(s) L∞ C.Time(s)

36 400 1.1419×10−12 0.270315 2.2197×10−12 0.434348
500 6.2560×10−12 1.212901 5.1609×10−12 1.193002
600 1.1793×10−12 1.119774 2.7161×10−12 0.904267
700 7.3037×10−12 1.541679 3.0467×10−12 1.403863

26 850 1.6612×10−11 1.114587 1.1432×10−11 1.473688
28 3.2163×10−12 1.106680 1.2373×10−11 1.363364
30 4.2749×10−12 1.425263 1.2061×10−11 1.471200
32 7.5047×10−12 2.004526 6.4597×10−12 1.563570
34 5.1494×10−12 1.505562 6.6363×10−12 1.666364

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.5

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.5

1
10

-12

(a)

0 0.2 0.4 0.6 0.8 1
10

-18

10
-16

10
-14

10
-12

10
-10

(b)
Figure 1: (a) Comparison of approximate and exact solutions in the subplot 1 and the LAbs

in the subplot 2 N = 700, MQ = 36 (Example 1). (b) Comparison of LAbs for different α
with N = 1000, MQ = 34 (Example 1).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 20 of 30

(a) (b)
Figure 2: (a) Graph of L∞ vs quadrature nodes MQ with N = 700 (Example 1). (b)
Graph of L∞ vs quadrature nodes N with MQ = 36 (Example 1).

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

10
-13

10
-12

10
-11

(a)

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

0.5

1

1.5

2

2.5

3

3.5

4

4.5
10

-12

(b)
Figure 3: (a) Graph of L∞ vs τ with MQ = 36 and N = 700 (Example 1). (b) Graph of
L∞ vs fractional order α with N = 700 and MQ = 36 (Example 1).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 21 of 30

(a)

Error Contour Plot

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2

-13.5

-13

-12.5

-12

-11.5

-11

(b)
Figure 4: (a) The graph shows L∞ error in τα plane with MQ = 36 and N = 700 (Example
1). (b) Contour plot of L∞ in τα plane with MQ = 36 and N = 700 (Example 1).

Example 2

In the second example, we consider the 2D version of (1)–(3) with λ1 = 1, λ2 = 0,
and exact solution u(x, τ) = (1− x2 − y2)τ.3. The performance of the proposed numerical
method is validated through detailed error analysis and computational tests. The L∞
error for various values of N and quadrature points MQ is shown in Table 2, showing
both high accuracy and computational efficiency. The numerical solution of example 2 is
shown in Figure 5a, and the surface plot presented in Figure 5b shows the absolute error
distribution demonstrating stable numerical performance.

Figure 6a demonstrates perfect quadrature convergence up to MQ = 36, after which the
error experiences a small rise. Spectral convergence with respect to spatial discretization
is evident in Figure 6b, where minor error oscillations appear at higher orders because
of round-off errors in Chebyshev differentiation matrices. Figure 7a shows the depen-
dence of L∞ on τ for 1.1 ≤ α ≤ 1.9, while Figure 7b shows its dependence on α for
τ = {0.1, 0.4, 0.7, 1}, both demonstrating consistently high accuracy. Figures 8a and 8b
present comprehensive error visualization using surface and contour plots in the ατ plane,
demonstrating the method’s robust stability across the entire parameter domain.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 22 of 30

Table 2: Errors norms for Example 2 with varying MQ, α, and N.

MQ N
α = 1.5 α = 1.75

L∞ C.Time(s) L∞ C.Time(s)

24 441 3.3135×10−9 0.525874 3.3135×10−9 0.430171
784 3.3137×10−9 1.371348 3.3135×10−9 1.370788
900 3.3135×10−9 1.790163 3.3138×10−9 2.238885
1089 3.3135×10−9 3.175910 3.3135×10−9 3.317322

26 1681 5.4329×10−10 9.138562 4.3794×10−10 9.125002
28 7.1942×10−10 9.830315 2.2237×10−10 9.696946
30 3.7555×10−10 10.442032 4.5945×10−10 10.396021
32 7.4227×10−10 11.273246 7.7937×10−10 11.177015
34 3.1544×10−10 13.683367 5.2071×10−10 12.233212

(a) (b)
Figure 5: (a) Comparison of approximate and exact solutions in the subplot 1 and the
LAbs in the subplot 2 N = 1681, MQ = 36 (Example 2). (b) Plot of LAbs of the method
N = 1681, MQ = 36 (Example 2).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 23 of 30

(a) (b)
Figure 6: (a) Graph of L∞ vs quadrature nodes MQ with N = 1681 (Example 2). (b)
Graph of L∞ vs quadrature nodes N with MQ = 36 (Example 2).

(a) (b)
Figure 7: (a) Graph of L∞ vs τ with MQ = 36 and N = 1681 (Example 2). (b) Graph of
L∞ vs fractional order α with N = 1681 and MQ = 36 (Example 2).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 24 of 30

(a)

-12

-11.5

-11

-10.5

-10

-9.5

(b)
Figure 8: (a) The graph shows the L∞ error in τα plane with MQ = 36 and N = 1681 (Ex-
ample 2). (b) Contour plot of L∞ in τα plane with MQ = 36 and N = 1681 (Example 2).

Example 3

In the third example, we consider the 3D version of (1)–(3) with λ1 = 1, λ2 = 0,
and exact solution u(x, τ) = exp(x + y + z)τ3. The L∞ error for various values of N and
quadrature points MQ is shown in Table 3, showing both high accuracy and computational
efficiency. The slice plots of numerical solution and absolute error computed usingN =
1331, MQ = 36, τ = 1, and α = 1.5 are presented in Figures 9a and 9b respectively. A
highly efficient performance is evident. The variation of L∞ vs MQ is presented in Figure
10a computed with N = 1728, α = 1.5, and τ = 1. Similarly, the variation of L∞ vs N is
shown in Figure 10b computed with MQ = 30, α = 1.5, and τ = 1. Figures 11a and 11b
presents the dependence of L∞ on τ and α respectively, both demonstrate high accuracy.
Further, Figures 12a show the error distribution in the αt plane. The contour slice plot
of absolute error is presented in Figure 12b. Overall, it is evident that the method has
the capability of solving fractional-order three-dimensional problems with high accuracy
without facing any time instability issues.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 25 of 30

Table 3: Errors norms for Example 2 with varying MQ, α, and N.

MQ N
α = 1.5 α = 1.75

L∞ C.Time(s) L∞ C.Time(s)

24 729 2.2140×10−9 1.265030 2.2140×10−9 1.013569
1331 2.2140×10−9 5.213977 2.2140×10−9 4.466417
2197 2.2140×10−9 17.880434 2.2140×10−9 15.620776
3375 2.2140×10−9 58.898524 2.2140×10−9 52.157271

26 1331 1.7283×10−10 5.531386 1.7283×10−10 4.843518
28 3.9919×10−11 6.134771 3.1216×10−11 5.339569
30 2.4134×10−11 6.425285 1.8025×10−11 5.591827
32 2.6654×10−11 5.969752 4.6171×10−11 6.005012
34 3.5469×10−11 6.343325 3.9908×10−11 6.308768

(a) (b)
Figure 9: (a) The slice plot of numerical solution with N = 1331, MQ = 36, τ = 1, and
α = 1.5 (Example 3). (b) The slice plot of LAbs with N = 1331, MQ = 36, τ = 1, and
α = 1.5 (Example 3).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 26 of 30

(a) (b)
Figure 10: (a) Graph of L∞ vs quadrature nodes MQ with N = 1728 (Example 3). (b)
Graph of L∞ vs quadrature nodes N with MQ = 30 (Example 3).

(a) (b)
Figure 11: (a) Graph of L∞ vs τ with MQ = 36 and N = 1728 (Example 3). (b) Graph
of L∞ vs fractional order α with N = 1728 and MQ = 36 (Example 3).



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 27 of 30

(a) (b)
Figure 12: (a) The graph shows the L∞ error in τα plane with MQ = 30 and N = 1728 (Ex-
ample 3). (b) Contour plot of L∞ in τα plane with MQ = 30 and N = 1728 (Example 3).

7. Conclusion

The paper develops the LT based CSCM method for numerical modelling of time-
fractional wave-diffusion equations including the MABC derivative. Unlike standard finite
difference methods, the proposed numerical method implements the LT and the numerical
inverse LT to efficiently handle the time-fractional derivative. It first utilizes the LT
to transform the considered problem into a time-independent inhomogeneous problem in
Laplace space. Then it employs the CSCM to discretize the spatial derivatives of the
transformed problem. Finally, it uses the improved Talbot method to recover the time-
domain solution.

Compared to conventional finite difference methods, the proposed LT-CSCM method
provides two key features: (i) elimination of computationally expensive convolution inte-
grals of fractional derivatives through LT; and (ii) unconditional stability independent of
time-stepping constraints. The features enable efficient and accurate long-time simulation
of diffusion-wave systems. The CSCM further improves the method’s efficiency for high-
dimensional problems, requiring fewer nodes while maintaining exponential convergence.
Numerical experiments confirm the LT based CSCM’s ability to handle multi-dimensional
diffusion-wave problems.

Looking forward, the robustness and efficiency of the LT based CSCM scheme make it
a strong candidate for simulating more complex fractional dynamical systems in applied
mathematics and engineering. Future work will focus on adapting this methodology to
solve fractional delay partial differential equations and coupled systems.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 28 of 30

Acknowledgements

The authors A. Aloqaily and N. Mlaiki would like to thank Prince Sultan University
for paying the publication fees for this work through TAS LAB.

Competing interests

There are no conflicting interests, according to the authors.

Author’s contributions

Each author contributed equally to the writing of this work, and they have all read
and approved the finished work.

Declarations

Ethical Approval Not applicable.
Funding This work did not receive any external funding.

References

[1] F. Mainardi. Fractional calculus: Theory and applications. Mathematics, 6(9):145,
2018.

[2] I. Podlubny. Fractional differential equations: an introduction to fractional deriva-
tives, fractional differential equations, to methods of their solution and some of their
applications, volume 198. Elsevier, 1998.

[3] Emad A Az-Zo’bi, Qais MM Alomari, Kallekh Afef, and Mustafa Inc. Dynamics of
generalized time-fractional viscous-capillarity compressible fluid model. Optical and
Quantum Electronics, 56(4):629, 2024.

[4] Mohammad A Al Zubi, Kallekh Afef, and Emad A Az-Zo’bi. Assorted spatial optical
dynamics of a generalized fractional quadruple nematic liquid crystal system in non-
local media. Symmetry, 16(6):778, 2024.

[5] S. G. Samko, A. A. Kilbas, and O. I. Marichev. Fractional integrals and derivatives:
Theory and applications. Gordon and Breach Science Publishers, 1993.

[6] F. Mainardi. Fractional calculus and waves in linear viscoelasticity. World Scientific,
2010.

[7] M. Al-Refai and D. Baleanu. On an extension of the operator with mittag-leffler
kernel. Fractals, 30(05):2240129, 2022.

[8] N. Shimizu and W. Zhang. Fractional calculus approach to dynamic problems of
viscoelastic materials. JSME International Journal Series C Mechanical Systems,
Machine Elements and Manufacturing, 42(4):825–837, 1999.



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 29 of 30

[9] K. Shah, H. Khalil, and R. A. Khan. Analytical solutions of fractional order diffusion
equations by natural transform method. Iranian Journal of Science and Technology,
Transactions A: Science, 42(3):1479–1490, 2018.

[10] M. Caputo. Linear models of dissipation whose q is almost frequency independent.
Annals of Geophysics, 19(4):383–393, 1966.

[11] M. A. Khan and A. Atangana. Modeling the dynamics of hepatitis e via the mab
derivative. Journal of Applied Mathematics and Computing, 55(1-2):345–358, 2017.

[12] F. Haq, K. Shah, G. ur Rahman, and M. Shahzad. Numerical solution of fractional
order smoking model via laplace adomian decomposition method. Alexandria Engi-
neering Journal, 57(2):1061–1069, 2018.

[13] K. Shah, M. A. Alqudah, F. Jarad, and T. Abdeljawad. Semi-analytical study of pine
wilt disease model with convex rate under caputo-febrizio fractional order derivative.
Chaos, Solitons & Fractals, 135:109754, 2020.

[14] F. Mainardi. The fundamental solutions for the fractional diffusion-wave equation.
Applied Mathematics Letters, 9(6):23–28, 1996.

[15] O. P. Agrawal. Solution for a fractional diffusion-wave equation defined in a bounded
domain. Nonlinear Dynamics, 29:145–155, 2002.

[16] J. Ren and Z. Z. Sun. Efficient numerical solution of the multi-term time fractional
diffusion wave equation. East Asian Journal on Applied Mathematics, 5(1):1–28, 2015.

[17] F. Liu, M. Meerschaert, R. McGough, P. Zhuang, and Q. Liu. Numerical methods for
solving the multi-term time-fractional wave-diffusion equation. Fractional Calculus
and Applied Analysis, 16(1):9–25, 2013.

[18] Y. N. Zhang, Z. Z. Sun, and X. Zhao. Compact alternating direction implicit scheme
for the two-dimensional fractional diffusion-wave equation. SIAM Journal on Numer-
ical Analysis, 50(3):1535–1555, 2012.

[19] R. Salehi. A meshless point collocation method for 2-d multi-term time fractional
diffusion-wave equation. Numerical Algorithms, 74:1145–1168, 2017.

[20] A. H. Bhrawy, E. H. Doha, D. Baleanu, and S. S. Ezz-Eldien. A spectral tau algorithm
based on jacobi operational matrix for numerical solution of time fractional diffusion-
wave equations. Journal of Computational Physics, 293:142–156, 2015.

[21] Y. Yang, Y. Chen, Y. Huang, and H. Wei. Spectral collocation method for the time-
fractional diffusion-wave equation and convergence analysis. Computers & Mathe-
matics with Applications, 73(6):1218–1232, 2017.

[22] M. H. Heydari, M. R. Hooshmandasl, F. M. Ghaini, and C. Cattani. Wavelets method
for the time fractional diffusion-wave equation. Physics Letters A, 379(3):71–76, 2015.

[23] F. A. Shah, Kamran, Z. A. Khan, F. Azmi, and N. Mlaiki. A hybrid collocation
method for the approximation of 2d time fractional diffusion-wave equation. AIMS
Mathematics, 9(10):27122–27149, 2024.

[24] M. Dehghan, M. Safarpoor, and M. Abbaszadeh. Two high-order numerical algo-
rithms for solving the multi-term time fractional diffusion-wave equations. Journal
of Computational and Applied Mathematics, 290:174–195, 2015.

[25] J. Y. Yang, J. F. Huang, D. M. Liang, and Y. F. Tang. Numerical solution of fractional
diffusion-wave equation based on fractional multistep method. Applied Mathematical



Kamran et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7095 30 of 30

Modelling, 38(14):3652–3661, 2014.
[26] B. L. Buzbee, G. H. Golub, and C. W. Nielson. On direct methods for solving

poisson’s equations. SIAM Journal on Numerical Analysis, 7(4):627–656, 1970.
[27] J. P. Boyd. Chebyshev and Fourier spectral methods. Dover Publications, 2001.
[28] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral methods: Fun-

damentals in single domains. Springer Science & Business Media, 2006.
[29] B. Dingfelder and J. A. C. Weideman. An improved talbot method for numerical

laplace transform inversion. Numerical Algorithms, 68(1):167–183, 2015.
[30] W. H. Huang, M. Samraiz, A. Mehmood, D. Baleanu, G. Rahman, and S. Naheed.

Modified atangana-baleanu fractional operators involving generalized mittag-leffler
function. Alexandria Engineering Journal, 75:639–648, 2023.

[31] R. Chawla, K. Deswal, D. Kumar, and D. Baleanu. A novel finite difference based
numerical approach for modified atangana-baleanu caputo derivative. AIMS Mathe-
matics, 7(9):17252–17268, 2022.

[32] P. Verma and M. Kumar. New existence, uniqueness results for multi-dimensional
multi-term caputo time-fractional mixed sub-diffusion and diffusion-wave equation on
convex domains. Journal of Applied Analysis and Computation, 11:1455–1480, 2021.

[33] L. N. Trefethen. Spectral methods in MATLAB. SIAM, Philadelphia, 2000.
[34] A. Shokri and S. Mirzaei. A pseudo-spectral based method for time-fractional

advection-diffusion equation. Computational Methods in Differential Equations,
8(3):454–467, 2020.

[35] B. D. Welfert. Generation of pseudospectral differentiation matrices I. SIAM Journal
on Numerical Analysis, 34(4):1640–1657, 1997.

[36] R. Baltensperger and M. R. Trummer. Spectral differencing with a twist. SIAM
Journal on Scientific Computing, 24(5):1465–1487, 2003.

[37] J. S. Green. The calculation of the time-responses of linear systems. PhD thesis,
Department of Applied Mathematics, Imperial College, London, 1955.

[38] A. Talbot. The accurate numerical inversion of laplace transforms. IMA Journal of
Applied Mathematics, 23(1):97–120, 1979.

[39] J. Weideman and L. N. Trefethen. Parabolic and hyperbolic contours for computing
the bromwich integral. Mathematics of Computation, 76(259):1341–1356, 2007.

[40] Steffen Börm, Lars Grasedyck, and Wolfgang Hackbusch. Introduction to hierarchical
matrices with applications. Engineering analysis with boundary elements, 27(5):405–
422, 2003.