EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6355 ISSN 1307-5543 – ejpam.com Published by New York Business Global Theoretical and Computational Analysis of Delay Volterra Integro-Differential Equations via Laplace Transform and Numerical Inversion Kamran1,∗, Nadeem Jan2, Muhammad Ishfaq Khan3, Ahmad Aloqaily4, Nabil Mlaiki4, Fady Hasan4 1 Department of Mathematics, Islamia College Peshawar, Peshawar 25120, Khyber Pakhtoonkhwa, Pakistan 2 Department of Natural Sciences and Humanities, University of Engineering & Technology Mardan, Mardan, Pakistan 3 College of Mechanics and Engineering Science, Hohai University, Nanjing 211100, China 4 Department of Mathematics and Sciences, Prince Sultan University, 11586 Riyadh, Saudi Arabia Abstract. Delay integro-differential equations (DIDEs) represent a significant class of integro- differential equations where state evolution depends on its past history. This paper presents a numerical approach for delay integro-differential equations (DIDEs), utilizing the Laplace trans- form (LT) and its inversion as the core methodology. Using the LT, the proposed technique begins by transforming the given equation into an algebraic equation in the Laplace domain. The resulting transformed equation is subsequently solved for the unknown function within the Laplace domain. Finally, two inversion methods, the Gauss-Hermite quadrature and the Weeks method, are used to invert the solution back to the time domain. Additionally, the existence and uniqueness of the solution are rigorously analyzed using functional analysis. To demonstrate the effectiveness of the methods, several examples from the literature are provided. The results obtained using the two techniques are compared and analyzed through tables and figures, highlighting their accuracy and computational efficiency. 2020 Mathematics Subject Classifications: 44A10, 65Rxx, 65R10 Key Words and Phrases: Delay integro-differential equation, Laplace transform, Existence, Uniqueness, Gauss-Hermite quadrature method, Weeks method. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6355 Email addresses: kamran.maths@icp.edu.pk (Kamran), maloqaily@psu.edu.sa (A. Aloqaily), nmlaiki@psu.edu.sa (N. Mlaiki), fhasan@psu.edu.sa (F. Hasan) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 2 of 22 1. Introduction Delay Volterra integro-differential equations (DVIDEs) represent an important class of functional equations that combines the complexities of neutral delay differential equations (NDDEs) with Volterra integro-differential (VID) systems. These equations are distin- guished by the integral terms that describe the cumulative effects of prior states and discrete and distributed delays in the dependent variables and their derivatives, so rep- resenting systems with memory and latency. DVIDEs have become important in recent research fields such as mathematical biology [1], engineering [2], physics [3], and economics [4], where dynamic systems exhibit temporal nonlocality and memory effects due to in- herent mathematical or physical constraints. For example, in biological models, DVIDEs are used to simulate epidemic spread with incubation times, population dynamics with gestational lags, and neural networks with synaptic delays. In engineering and physical sciences, they are used to model anomalous diffusion, viscoelastic materials, and time- lag feedback control systems. In economics, these equations assist in analyzing complex systems with delayed feedback loops, including market reactions, supply chains, and in- vestment decisions, establishing a strong basis to analyze non-instantaneous causality in changing situations. The numerical analysis of solutions to DVIDEs is a progressing topic of research. For example, the authors of [5] applied the multi-step Adams-Moulton method to solve the neutral DVIDEs. Mansouri and Azimzadeh [6] used the Bernstein polynomial approach to solve the DVIDEs by transforming the delay terms and the integral components into a polynomial framework. In [7], the authors have explored the explicit and implicit Runge- Kutta methods for solving the neutral DVIDEs. They analyzed the convergence properties of the scheme required at each step as well as the overall convergence of the numerical scheme. Zaidan [8] presented a method for approximating the solution of linear VIDEs. The suggested approach employs Galerkin’s approach, with Bernstein polynomials as the basis functions. The authors of [9] explored different aspects of delay Volterra functional integral equations, concentrating on the existence, uniqueness, and regularity of solutions. Brunner [10] presented different results on the analysis of local and global super conver- gence orders in collocation methods for delay Volterra functional differential equations. Bellour and Bousselsal [11] studied the numerical solution of DVIDEs using the Taylor collocation method. The authors of [12] investigated the error behaviour of linear multi- step methods applied to singularly perturbed DVIDEs. Amirali and Acar [13] introduced a novel numerical technique for solving neutral DVIDEs. Further, their numerical scheme obtained a second-order convergence. Rihan et al. [14] presented a new numerical scheme for solving the DVIDEs; the authors used the implicit Runge-Kutta method for approxi- mating the differential operator and Boole’s quadrature rule for the integral operator. Even though traditional numerical methods offer viable solutions to DVIDEs, they often struggle to handle complex integrals and decay terms effectively. In this context, the Laplace transform (LT) emerges as a promising method, providing an elegant approach to simplify these equations by converting them into the Laplace domain, where algebraic manipulation becomes easy. The use of LT facilitates the analysis of systems with delays, Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 3 of 22 hereditary properties, or distributed memory effects, as these intricate processes become more feasible in transformed space. Moreover, the LT method is particularly effective in solving initial value problems because it integrates the initial conditions directly into the transformed equation, eliminating the need for additional steps. Multiple studies have shown the efficiency of the LT method in solving equations related to population dynamics, viscoelasticity problems, control theory, etc ([15–18]). However, the success of this method is largely dependent on the availability of effective numerical inversion methods to extract the time-domain solution from the Laplace domain solution [19, 20]. In the literature, many methods have been developed for numerical inversion of LT. However, in this article, we use the Gauss-Hermite quadrature (GHQ) method [21] and the Weeks (WK) method [22]. The structure of this paper is organized as follows: Section 2 presents the fundamen- tal definitions. Section 3 explores the existence and uniqueness results of the solution. Section 4 provides a detailed description of the methods, followed by Section 5, which includes numerical examples that illustrate the efficiency of the methods. Finally, Section 6 concludes the article by presenting key insights and outlining potential future directions. 2. Preliminaries Let T = [0, 1]. For any function u ∈ C(T , R), the norm ∥u∥∞ is defined as ∥u∥∞ = sup t∈T |u(t)|. Definition 1. Let u(t) be a real-valued function that is piecewise continuous for t > 0, and assume that it is of exponential order. Under these conditions, the LT of u(t) exists and is given as: û(z) = L {u(t)} = ∫ ∞ 0 e−ztu(t)dt. Theorem 1. [3] Let ϕ(t),ϕ′(t) be continuous on the interval [−σ, 0]; then the LT of u(t−σ) and u′(t− σ) are given as: L {u(t− σ)} = ϕ̄(z) + e(−zσ)L {u(t)}, and L {u′(t− σ)} = ¯̄ϕ(z) + e(−zσ)L {u′(t)}, where ϕ̄(z) = ∫ 0 −σ e−z(t+σ)ϕ(t)dt, ¯̄ϕ(z) = ∫ 0 −σ e−z(t+σ)ϕ′(t)dt. 3. Existence and uniqueness of solution The analysis of the existence and uniqueness of the solution is crucial for ensuring that the problem at hand is well-posed. By establishing these characteristics, the model Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 4 of 22 becomes reliable and consistent, ensuring that a solution exists under specific conditions and is unique. This section offers a rigorous framework to show that the problem con- sidered admits a unique solution, setting a strong foundation for subsequent numerical investigations. In this study, the following DVIDE is analyzed: λ1 du(t) dt +λ2 du(t− σ) dt = λ3f(t)+λ4u(t)+λ5u(t−σ)+ ∫ t 0 K(t, x) (λ6u(x) + λ7u(x− σ)) dx, t ∈ T , (1) u(t) = ϕ(t), t ∈ [−σ, 0], u(0) = u0, where λ1, λ2, λ3, λ4, λ5, λ6, λ7 are constants, the kernel K(t, x) = K(t−x) is a smooth function, f(t) is a given function, ϕ is a delay condition and u(t) is the unknown function to be determined. Applying integration techniques to the problem in Eq (1) yields the following solution. u(t) = u0 + 1 λ1 [ λ2(u(−σ) − u(t− σ)) + ∫ t 0 λ3f(τ) dτ + ∫ t 0 λ4u(τ) dτ + ∫ t 0 λ5u(τ − σ) dτ + ∫ t 0 ∫ τ 0 K(x, τ) (λ6u(x) + λ7u(x− σ)) dx dτ ] . Now, we consider the operator T : C([0, 1]) → C([0, 1]), defined by: (Tu)(t) = u0 + λ2 λ1 u(−σ) + λ2 λ1 u(t− σ) + λ3 λ1 ∫ t 0 f(s)ds+ λ4 λ1 ∫ t 0 u(s)ds+ λ5 λ1 ∫ t 0 u(s− σ)ds + λ6 λ1 ∫ t 0 (∫ s 0 K(x, s)u(x)dx ) ds+ λ7 λ1 ∫ t 0 (∫ s 0 K(x, s)u(x− σ)dx ) ds, for all t ∈ [0, 1] and u ∈ C([0, 1]), where u(t) = ϕ(t), for t ∈ [−σ, 0]. Lemma 1. A function u ∈ C([0, 1]) is a solution to the integral equation u = Tu if and only if u ∈ C1([0, 1]) and satisfies the DVIDE in Eq.(1) with the initial condition u(t) = ϕ(t) for t ∈ [−σ, 0]. Proof. Suppose u ∈ C([0, 1]) satisfies u = Tu. Substituting u(t) = ϕ(t) for t ∈ [−σ, 0], the integral equation matches the form obtained by integrating Eq.(1). If u ∈ C1([0, 1]), differentiate both sides of the integral equation with respect to t. Using the fundamental theorem of calculus and the continuity of K, we recover Eq.(1): λ1 du(t) dt +λ2 du(t− σ) dt = λ3f(t)+λ4u(t)+λ5u(t−σ)+ ∫ t 0 K(t, x) (λ6u(x) + λ7u(x− σ)) dx. In contrast, if u ∈ C1([0, 1]) satisfies Eq.(1), integrating both sides from 0 to t and applying the initial condition u(0) = u0 yields the integral equation. Thus, the solutions C1([0, 1]) in the DVIDE correspond to the solutions in C([0, 1]) to the integral equation that are differentiable. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 5 of 22 Theorem 2. The problem defined in Eq (1) has at least one solution. Proof. The proof involves several steps. Step 1: Continuity of the Operator T For any continuous function u, Tu is also continu- ous. Indeed, consider un, u ∈ C([0, 1]). We have |Tun(t) − Tu(t)| = ∣∣∣∣∣λ2λ1 (un(−σ) − u(−σ)) + λ2 λ1 (un(t− σ) − u(t− σ)) + λ4 λ1 ∫ t 0 ( un(s) − u(s) ) ds + λ5 λ1 ∫ t 0 ( un(s− σ) − u(s− σ) ) ds+ λ6 λ1 ∫ t 0 (∫ s 0 K(x, s) ( un(x) − u(x) ) dx ) ds + λ7 λ1 ∫ t 0 (∫ s 0 K(x, s) ( un(x− σ) − u(x− σ) ) dx ) ds ∣∣∣∣∣ ≤ λ2 λ1 |un(−σ) − u(−σ)| + λ2 λ1 |un(t− σ) − u(t− σ)| + λ4 λ1 ∫ t 0 |un(s) − u(s)|ds + λ5 λ1 ∫ t 0 |un(s− σ) − u(s− σ)|ds+ λ6 λ1 ∫ t 0 (∫ s 0 (|K(x, s)||un(x) − u(x)|) dx ) ds + λ7 λ1 ∫ t 0 (∫ s 0 (|K(x, s)||un(x− σ) − u(x− σ)|) dx ) ds, taking the supremum norm on both sides yields: ∥Tun − Tu∥∞ ≤ 2λ2 λ1 ∥un − u∥∞ + λ4 λ1 ∫ t 0 ∥un − u∥∞ds+ λ5 λ1 ∫ t 0 ∥un − u∥∞ds + λ6 λ1 ∫ t 0 (∫ s 0 (N∥un − u∥∞) dx ) ds+ λ7 λ1 ∫ t 0 (∫ s 0 (N∥un − u∥∞) dx ) ds, where N = max |K(x, s)|. Simplifying gives ∥Tun − Tu∥∞ ≤ ( 2λ2 λ1 + λ4 λ1 + λ5 λ1 + λ6N 2λ1 + λ7N 2λ1 ) ∥un − u∥∞. Since u is continuous, it follows that ∥Tun − Tu∥∞ → 0 as n→ ∞. Thus, the operator T is continuous. Step 2: Boundedness of T We show that the operator maps bounded sets into bounded sets. For u ∈ Bη1 = {u ∈ C([0, 1]) : ∥u∥∞ ≤ η1}, we have |Tu(t)| ≤ |u0| + λ2 λ1 |u(−σ)| + λ2 λ1 |u(t− σ)| + λ3 λ1 ∫ t 0 |f(s)|ds+ λ4 λ1 ∫ t 0 |u(s)|ds+ λ5 λ1 ∫ t 0 |u(s− σ)|ds + λ6 λ1 ∫ t 0 (∫ s 0 |K(x, s)||u(x)|dx ) ds+ λ7 λ1 ∫ t 0 (∫ s 0 |K(x, s)||u(x− σ)|dx ) ds. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 6 of 22 Using boundedness of f and u, we obtain ∥Tu∥∞ ≤ |u0| + 2λ2 λ1 η1 + λ3 λ1 M + λ4 + λ5 λ1 η1 + (λ6 + λ7)Nη1 2λ1 =: η2, where M = max |f(s)| and N = max |K(x, s)|. This shows that T maps bounded sets into bounded sets. Step 3: Equicontinuity of T To show equicontinuity, for t, t0 ∈ [0, 1], we have |Tu(t) − Tu(t0)| = ∣∣∣∣(u0 + λ2 λ1 u(−σ) + λ2 λ1 u(t− σ) + λ3 λ1 ∫ t 0 f(s)ds+ λ4 λ1 ∫ t 0 u(s)ds+ λ5 λ1 ∫ t 0 u(s− σ)ds + λ6 λ1 ∫ t 0 (∫ s 0 K(x, s)u(x)dx ) ds+ λ7 λ1 ∫ t 0 (∫ s 0 K(x, s)u(x− σ)dx ) ds ) − ( u0 + λ2 λ1 u(−σ) + λ2 λ1 u(t0 − σ) + λ3 λ1 ∫ t0 0 f(s)ds+ λ4 λ1 ∫ t0 0 u(s)ds+ λ5 λ1 ∫ t0 0 u(s− σ)ds + λ6 λ1 ∫ t0 0 (∫ s 0 K(x, s)u(x)dx ) ds+ λ7 λ1 ∫ t0 0 (∫ s 0 K(x, s)u(x− σ)dx ) ds )∣∣∣∣ ≤ λ2 λ1 |u(t− σ) − u(t0 − σ)| + λ3 λ1 ∫ t t0 |f(s)|ds+ λ4 λ1 ∫ t t0 |u(s)|ds+ λ5 λ1 ∫ t t0 |u(s− σ)|ds + λ6 λ1 ∫ t t0 (∫ s 0 |K(x, s)||u(x)|dx ) ds+ λ7 λ1 ∫ t t0 (∫ s 0 |K(x, s)||u(x− σ)|dx ) ds ≤ Rλ2 λ1 |t− t0| + λ3 λ1 ∫ t t0 |f(s)|ds+ λ4 λ1 ∫ t t0 |u(s)|ds+ λ5 λ1 ∫ t t0 |u(s− σ)|ds + λ6N λ1 ∫ t t0 (∫ s 0 |u(x)|dx ) ds+ λ7N λ1 ∫ t t0 (∫ s 0 |u(x− σ)|dx ) ds. Simplifying, ∥Tu(t)−Tu(t0)∥∞ ≤ Rλ2 λ1 ∥t−t0∥∞+ λ3 λ1 M(t−t0)+ λ4 + λ5 λ1 η1(t−t0)+ (λ6 + λ7)Nη1 2λ1 (t2−t20), where |u(t−σ)−u(t0−σ)| < R|t− t0|. Hence ∥Tu(t)−Tu(t0)∥∞ → 0 as t→ t0. Therefore by Arzelà-Ascoli Theorem [23], the operator T is completely continuous. Step 4: A Priori Bound Define the set ω = {u ∈ C([0, 1]) : u = ϵTu, 0 < ϵ < 1}. For u ∈ ω, ∥u∥∞ ≤ ϵ∥Tu∥∞ ≤ ϵη2. Which demonstrates the boundedness of ω. Utilizing Schaefer’s fixed point theorem, it follows that the operator T possesses at least one fixed point. As a result, the problem described in Eq.(1) guarantees the existence of at least one solution. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 7 of 22 Theorem 3. The problem in Eq.(1) has a unique solution in C([0, 1]). Proof. For u1, u2 ∈ C([0, 1]), compute: |Tu1(t)−Tu2(t)| ≤ λ2 λ1 |u1(−σ)−u2(−σ)|+λ2 λ1 |u1(t−σ)−u2(t−σ)|+λ4 λ1 ∫ t 0 |u1(s)−u2(s)| ds + λ5 λ1 ∫ t 0 |u1(s− σ) − u2(s− σ)| ds+ λ6 λ1 ∫ t 0 ∫ s 0 |K(x, s)||u1(x) − u2(x)| dx ds + λ7 λ1 ∫ t 0 ∫ s 0 |K(x, s)||u1(x− σ) − u2(x− σ)| dx ds. Since u1(−σ) = u2(−σ) = ϕ(−σ), the first term is zero. When t − σ ≥ 0, |u1(t − σ) − u2(t− σ)| ≤ ∥u1 − u2∥∞; similarly for s− σ ≥ 0, x− σ ≥ 0. Thus: ∥Tu1 − Tu2∥∞ ≤ ( 2λ2 λ1 + λ4 + λ5 λ1 + (λ6 + λ7)N 2λ1 ) ∥u1 − u2∥∞. This constant may not be less than 1, so T may not be a contraction. For the Volterra-type equation, consider T k. For k ≥ ⌈ 1 σ ⌉, delay terms like u(t − kσ) in T ku1 − T ku2 satisfy t − kσ ≤ 1 − kσ ≤ 0, since t ≤ 1, so they equal ϕ(t − kσ), making their differences zero. Therefore: ∫ t 0 ∫ s1 0 · · · ∫ sk−1 0 |u1(x) − u2(x)| dx dsk−1 · · · ds1 ≤ 1 k! ∥u1 − u2∥∞. Thus: ∥T ku1 − T ku2∥∞ ≤ ( λ4 + λ5 + (λ6+λ7)N 2 λ1 1 k! ) ∥u1 − u2∥∞. Choose k such that: λ4 + λ5 + (λ6+λ7)N 2 λ1 1 k! < 1, so T k is a contraction. By the Banach fixed-point theorem, T has a unique fixed point in C([0, 1]). By Lemma 1, the DVIDE has a unique solution in C1([0, 1]). 4. Methodology This section outlines the numerical scheme devised for solving DIDEs. The key steps involved are: (i) applying the LT to transform the DIDE to an algebraic equation in the LT domain; (ii) the resulting algebraic equation is solved in the transformed domain; and (iii) utilizing the inverse LT to retrieve the solution in the time domain. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 8 of 22 4.1. Laplace transform application By applying the LT to Eq. (1), we obtain: L { λ1 du(t) dt + λ2 du(t− σ) dt = λ3f(t) + λ4u(t) + λ5u(t− σ) + ∫ t 0 K(t, x) (λ6u(x) + λ7u(x− σ)) dx } , or L { λ1 du(t) dt } + L { λ2 du(t− σ) dt } = L {λ3f(t)} + L {λ4u(t)} + L {λ5u(t− σ)} + L {∫ t 0 K(t, x) (λ6u(x) + λ7u(x− σ)) dx } , which leads to λ1{zû(z) − u(0)} + λ2{ ¯̄ϕ(z) + e−zσ {zû(z) − u(0)} = λ3f̂(z) + λ4û(z) + λ5{ϕ̄(z) + e−zσû(z)} + K̂(z) { λ6û(z) + λ7{ϕ̄(z) + e−zσû(z)} } , which implies{ λ1z + λ2ze −zσ − λ4 − λ5e −zσ − K̂(z)(λ6 + λ7e −zσ) } û(z) = λ1u0 − λ2 ¯̄ϕ(z) + λ2e −zσu0 + λ3f̂(z) + λ5ϕ̄(z) + λ7K̂(z)ϕ̄(z), leads to û(z) = λ1u0 − λ2 ¯̄ϕ(z) + λ2e −zσu0 + λ3f̂(z) + λ5ϕ̄(z) + λ7K̂(z)ϕ̄(z) λ1z + λ2ze−zσ − λ4 − λ5e−zσ − K̂(z)(λ6 + λ7e−zσ) , hance, the solution of problem in Eq. (1) can be obtained by applying the inverse LT as: u(t) = 1 2πi ∫ γ+i∞ γ−i∞ eztû(z)dz γ > γ0. (2) The function û(z) is presumed to be analytic within the half-plane defined by ℜ(z) > γ0, where γ0 represents the convergence abscissa. Our objective is to evaluate the original solution u(t) for one or more values of t > 0. Typically, a contour C is selected connecting the points γ− i∞ and γ+ i∞ along a line ℜ(z) = γ. Nevertheless, using Cauchy’s theorem allows for the deformation of C into a shape that is suitable for the approximation of (2). Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 9 of 22 4.2. Gauss-Hermite-quadrature method The GHQ method is a numerical integration technique that approximates the integral of functions multiplied by a Gaussian weight function exp(−κ2). It is especially useful for integrals of the form ∫ ∞ −∞ exp(−t2)g(t)dt (3) The GHQ method approximates Eq. (3) by evaluating the function u(t) at certain nodes and weighting them appropriately[21]∫ ∞ −∞ exp(−t2)g(t)dt ≈ ν∑ ϱ=1 ηϱg(tϱ), (4) where • tϱ represents the roots of Hermite polynomial Hν(t), of degree ν, • ηϱ are the associated weights, The ηϱ are uniquely determined by the fact that, for a polynomial Hν(t) of degree at most ν − 1, the approximation becomes exact. GHQ quickly converges for smooth Hν(t). The following theorem in [24, 25] clarifies this point: Theorem 4. The error in (4) can be written as follows: Rν(g) = 1 2πi ∫ C ψν(z)g(z)dz, (5) where ψν(z) = Sν(z) Hν(z) , Sν(z) = ∫ ∞ −∞ e−t2 Hν(t) z − t dt. (6) The contour C in (5) starts at z = −∞, and loops around the zeros of Hν , and terminates at z = ∞, assuming f(z) is analytic within the region enclosed by C. The Hermite function of the second kind, Sν(z), defined in (6), exhibits rapid decay as |z| → ∞, a feature that facilitates the convergence of integral expressions involving it. Conversely, Hν(z) grows exponentially as |z| → ∞, which constraints the function ψν(z) in (5) to remain bounded, provided the contour C is carefully chosen to account for the growth and singularities of f(z). The accuracy and stability of the method depend on the growth of Hν(z) and decay of Sν(z). Furthermore, related studies (e.g.,) highlight the importance of contour deformation to minimize numerical errors and indicate that high accuracy is achieved when singular- ities of f(z) are sufficiently distant from the real axis. This characteristic is especially advantageous for guaranteeing accuracy in situations where singularities might affect the contour’s path. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 10 of 22 By selecting t = ςκ, with ς ∈ R, allows for a simple integration process while main- taining high accuracy, as suggested in past analysis of similar integral representations in orthogonal polynomial theory:∫ ∞ −∞ e−t2g(t)dt = ς ∫ ∞ −∞ e−ς2κ2 g(ςκ)dκ = ς ∫ ∞ −∞ e−κ2 eκ 2 e−ς2κ2 g(ςκ)dκ. (7) Now, we can utilize the rule (4) to the function F(κ) = eκ 2(1−ς2)f(ςκ). As a result the function F̂(κ), will exhibit singularities that are fartherfrom the x-axiscompared to those of f(z). This increased separation of singularities has the capability to improve the numerical accuracy. To evaluate (2), we consider the following contour: C: z = ϑ(1 + iγ)2, where γ ∈ (−∞,∞), and ϑ > 0. (8) The contour described in (8) intersects the x − axis at z = ϑ and the y − axis at z = ±2ϑi. Utilizing (8) in (2), we obtain∫ C û(z)ezκdz = ∫ ∞ −∞ ez(γ)κû(z(γ))z′(γ)dγ, the technique in (7), can be applied using γ = ςσ as: u(κ) = ς 2πi ∫ ∞ −∞ e−σ2 eσ 2+z(ςσ)κû(z(ςσ))z′(ςσ)dσ = ∫ ∞ −∞ e−σ2 g(σ)dσ, leads to u(κ) = ∫ ∞ −∞ e−σ2 g(σ)dσ, (9) where g(σ) = ς 2πi eσ 2+z(ςσ)κû(z(ςσ))z′(ςσ). Applying the GHQ to (9) yields uApp(κ) ≈ 2Re  p∑ ϱ=1 ηϱg(σϱ)  . Here, σϱ represents the positive roots of Hν , and for even ν, p = ν/2, while for odd ν, p = (ν + 1)/2. 4.2.1. Error analysis The error analysis in the GHQ rule can be analyzed using the following theorem: Theorem 5. [26] Let g(σϱ) be the transformed integrand in the GHQ rule. If g(σϱ) is 2ν-times continuously differentiable on (−∞,∞), then the error in Eν in the quadrature approximation satisfies: |Eν | ≤ C 2ν! ∥g(2ν)∥∞, where C is a constant depending on ν. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 11 of 22 Proof. The GHQ rule is exact for polynomials of degree up to 2ν − 1. We begin by expanding g(σϱ) using a Taylor series as: g(σϱ) = P2ν−1(σϱ) +R2ν(σϱ), where R2ν(σϱ) = g(2ν)(ξ) (2ν)! (σϱ − σϱ0)2ν . The error Eν arises from the remainder term: Eν = I(R2ν) −Qν(R2ν). Further, we have |R2ν(σϱ)| ≤ ∥g(2ν)∥∞ (2ν)! |σϱ − σϱ0 |2ν . Now, integrating the remainder term we have |Eν | ≤ ∥g(2ν)∥∞ (2ν)! ∫ ∞ −∞ e−σ2 ϱ |σϱ − σϱ0 |2νdσϱ. Let C = ∫∞ −∞ e−σ2 ϱ |σϱ − σϱ0 |2νdσϱ. Then, we obtain |Eν | ≤ C (2ν)! ∥g(2ν)∥∞. 4.3. Weeks method The Weeks method is acknowledged as one of the most straightforward, accurate, and reliable methods for the numerical inversion of the LT, as long as the parameters involved in the Laguerre series expansion are properly chosen. Unlike the trapezoidal rule and Talbot’s technique, the Weeks method offers a significant advantage through its function expansion, specifically using the Laguerre series. This expansion allows the determination of unknown coefficients for any given û(z) ensuring high accuracy. A key feature of Weeks’ approach is its ability to include a complex parameter z defined as z = χ+ iϑ, ϑ ∈ R. In addition to improving accuracy, this parametrization offers more flexibility to a wide range of functions with different singularity structures. Studies like those carried out by Abate and Whitt [27] have shown the efficiency of the Weeks method in minimizing the oscillations during the inversion technique, especially when compared to other inversion methods. It is a preferred choice in a wide range of applied and theoret- ical contexts due to its versatility, which is demonstrated by its ability to adapt to both smooth and non-smooth functions. Thus, by accurately tuning the parameters χ and ξ, the Weeks technique guarantees a solid and efficient technique, leading to the formulation as follows: Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 12 of 22 un(t) = eχt 2π ∫ ∞ −∞ eitξû(χ+ iξ)dξ. (10) The transform function û(χ+ iξ) is expanded using the Laguerre series of the form: û(χ+ iξ) = ∞∑ ȷ=−∞ αȷ (−ϖ + iξ)ȷ (ϖ + iξ)ȷ+1 , ϖ > 0, ξ ∈ R. (11) Using (11) in (10), resulting in: un(t) = eχt 2π ∞∑ ȷ=−∞ αȷδȷ(t;ϖ), (12) where δȷ(t;ϖ) = ∫ ∞ −∞ eitξ (−ϖ + iξ)ȷ (ϖ + iξ)ȷ+1 dξ. (13) For evaluating the Fourier integral, residues offer a potent tool, particularly when dealing with complex-valued functions. For t > 0, the evaluation yields: δȷ(t;ϖ) = { 2πe−ϖtLȷ(2ϖt), ȷ ≥ 0, 0, ȷ < 0. (14) The ȷth degree Laguerre polynomials, Lȷ(t), play an important role in expansions like those used in the Weeks method. The Lȷ(t) are defined as: Lȷ(t) = et ȷ! dȷ dtȷ (e−ttȷ). (15) In this contexts, χ > χ0, where χ0 denotes the abscissa of convergence. The parameters χ,ϖ represent positive real numbers, while αȷ are the unknown coefficients that must be determined. C(φ) = 2ϖ 1 − φ û ( χ+ 2ϖ 1 − φ −ϖ ) = ∞∑ ȷ=0 αȷφ ȷ, |φ| < R, (16) where R represents the radius of convergence of the Maclaurin series in Eq. (16). The coefficients αȷ are calculated as follows: αȷ = 1 2πi ∫ |φ|=1 C(φ) φȷ+1 dφ = 1 2π ∫ π −π C(eiζ)e−iȷζdζ. (17) The integral in Eq. (17) is equivalent to Cauchy’s integral formula, which is crucial for evaluating contour integrals. The integral can be numerically approximated as follows: α̃ȷ = e−iȷλ/2 2n n−1∑ j=−n C(eiζj+1/2)e−iȷζj , ȷ = 0, 1, 2, ..., n− 1, (18) where ζj = jλ, λ = π n . Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 13 of 22 4.3.1. Error analysis Error analysis plays an important role in assessing the reliability and precision of numeri- cal methods, especially for algorithms designed to perform numerical inversion of the LT. Analyzing the cause and propagation of errors allows the assessment of the efficiency of selected methods and the identification of potential constraints. The Weeks method, a widely used method for numerical LT inversion, adds several error factors that can affect the accuracy of the results. These errors result from the numerical computation of the coefficient, the truncation of infinite series, and the challenging process of numerical inver- sion of the LT. By comprehending and quantifying these errors, researchers can enhance their approach and ensure the accuracy of their solutions. Weideman [22] investigated the error in the Weeks method and provided detailed observations. u(t) = exp(χt) ∞∑ ȷ=0 αȷexp(−ϖt)Lȷ(2ϖt). (19) Three primary sources of error were identified. • 1st : The first source arise from truncation of the Laguerre series, • 2nd : The second source stems from computing the unknown Laguerre coefficients numerically, • 3rd : The third source of error is the numerical inversion of the LT. The solution incorporating these errors is expressed as follows: ũ(t) = exp(χt) n−1∑ ȷ=0 α̃ȷ(1 + δȷ)exp(−ϖt)Lȷ(2ϖt), (20) where δȷ represents the relative error in the floating-point representation of the coefficients, i.e., fl(α̃ȷ) = α̃ȷ(1 + δȷ). Subtracting (20) from (19) and letting ∑∞ ȷ=0 |αȷ| <∞, the total error can be bounded as: |u(t) − ũ(t)| ≤ e(χt) ( erT + erD + erC ) , where erT is the truncation error bound defined as: erT = ∞∑ ȷ=n |αȷ|, erD is the discretization error bound defined as: erD = n−1∑ ȷ=0 |αȷ − α̃ȷ|, Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 14 of 22 and erC is the conditioning error bound defined as: erC = χ n−1∑ ȷ=0 |α̃ȷ|. The machine round-off unit, denoted by δ satisfies max0≤ȷ≤n−1|δȷ| ≤ δ.. In partic- ular, we assume |exp(−ϖt)Lȷ(2ϖt)| ≤ 1. To simplify the error analysis, we neglect the discretization error erD and focus on erT and erC [22]. The bounds provided in [22] are: erT ≤ ȷ(t) tn(t− 1) , erC ≤ δ tȷ(t) (t− 1) , where1 < t < R. Therefore, the overall error estimate is given by: errorest ≤ ȷ(t) tn(t− 1) + δ tȷ(t) t− 1 . (21) 5. Application This section provides the computational results for four different problems solved using the proposed numerical schemes. Extensive research has addressed these problems, and their numerical solutions are well documented, providing a solid basis for comparative analysis. To assess the accuracy of the numerical schemes, two error matrices are used: the Lin and the Lrms, defined as follows: Lin = max 1≤j≤n |u(tj) − un(tj))|, and Lrms = √∑n j=1(u(tj) − un(tj))2 n , where u(t) represent the exact solution and un(t) represent the numerical solution of the considered problem. 5.1. Problem 1 In the first problem, we consider Eq.(1), with parameters λ1 = λ4 = λ5 = λ7 = 0, λ2 = λ3 = λ6 = 1, σ = 1, kernel function K(t, x) = t− x, and source term f(t) = 1 − t3 6 , with ϕ(t) = t and exact solution u(t) = t. The Lin and Lrms are provided in Table 1, demonstrating the convergence of both schemes. It is observed that WM achieves high accuracy as compared to the GHQ. However, when comparing computational times, the GHQ demonstrates better efficiency. The numerical results of both methods are also compared with the results from [28] and [29], showing that the proposed numerical schemes provide better accuracy. The exact and numerical solutions for this problem are illustrated in Fig. 1a, while the Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 15 of 22 comparison of Lin and Lrms for the two numerical schemes is shown in Fig. 1b. These comparisons clearly indicate that the WM method outperforms the GHQ in terms of accuracy. Table 1: The Lin and Lrms corresponding to Problem 1 using GHM and WM 1.3 GHM WM t Lin Lrms Lin Lrms 0.1 3.61×10−14 8.08×10−15 2.78×10−17 3.93×10−18 0.2 7.21×10−14 1.61×10−14 8.33×10−17 9.21×10−18 0.3 1.09×10−13 2.45×10−14 1.11×10−16 1.44×10−17 0.4 1.44×10−13 3.22×10−14 1.11×10−16 1.82×10−17 0.5 1.82×10−13 4.06×10−14 1.11×10−16 1.90×10−17 0.6 2.16×10−13 4.83×10−14 1.11×10−16 2.20×10−17 0.7 2.56×10−13 5.72×10−14 2.22×10−16 3.13×10−17 0.8 2.89×10−13 6.46×10−14 2.22×10−16 3.84×10−17 0.9 3.27×10−13 7.31×10−14 2.22×10−16 3.99×10−17 1 3.62×10−13 8.10×10−14 6.66×10−16 7.77×10−17 CPU(s) 0.001242 1.121635 [28] 1.88×10−8 1.571 [29] 9.28×10−8 —- 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Exactsol numericsol (a) 0 0.2 0.4 0.6 0.8 1 10-18 10-17 10-16 10-15 10-14 10-13 10-12 L in (GHM) L in (WM) L rms (GHM) L rms (WM) (b) Figure 1: (a) Comparison of approximate and exact solutions for problem 1. (b) The variation of error norms vs t for problem 1. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 16 of 22 5.2. Problem 2 In the second problem, we consider Eq.(1), with parameters λ2 = λ5 = λ6 = 0, λ1 = λ3 = λ4 = λ7 = 1, σ = 1, the kernel function K(t, x) = 1, ϕ(t) = et, and the source term f(t) = e−1(1 − et). The exact solution for this problem is u(t) = et. The problem is solved using the two numerical schemes, the GHQ and WM. The Lin and Lrms for both methods are presented in Table 2, highlighting their convergence. It is observed that WM provides high accuracy compared to the GHQ. However, in terms of computational efficiency, the GHQ demonstrates better performance. Additionally, the numerical results of both methods are compared with those reported in [30], showing that the proposed numerical schemes yield superior accuracy. The exact and numerical solution for this problem is illustrated in Fig. 2a, while the comparison of Lin and Lrms for the two inversion schemes is shown in Fig. 2b. These results clearly indicate that the WM outperforms the GHQ in terms of accuracy. Table 2: The Lin and Lrms corresponding to Problem 2 using the GHM and WM. 1.3 GHM WM t Lin Lrms Lin Lrms [30] 0.1 1.95×10−13 4.35×10−14 2.22×10−16 2.22×10−17 1.98×10−5 0.2 1.49×10−12 3.33×10−13 2.22×10−16 3.14×10−17 2.15×10−5 0.3 7.15×10−12 1.60×10−12 2.22×10−16 3.14×10−17 2.34×10−5 0.4 2.71×10−11 6.06×10−12 4.44×10−16 5.44×10−17 2.56×10−5 0.5 8.85×10−11 1.98×10−11 4.44×10−16 5.87×10−17 2.81×10−5 0.6 2.59×10−10 5.80×10−11 4.44×10−16 5.87×10−17 3.08×10−5 0.7 7.01×10−10 1.57×10−10 4.44×10−16 5.87×10−17 3.40×10−5 0.8 1.78×10−09 3.98×10−10 8.88×10−16 1.06×10−16 3.75×10−5 0.9 4.28×10−09 9.57×10−10 8.88×10−16 1.39×10−16 4.16×10−5 1 9.86×10−09 2.20×10−09 8.88×10−16 1.65×10−16 4.61×10−5 CPU(s) 0.167872 1.181213 0.0523884 Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 17 of 22 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 Exactsol numericsol (a) 0 0.2 0.4 0.6 0.8 1 10-18 10-16 10-14 10-12 10-10 10-8 L in (GHM) L in (WM) L rms (GHM) L rms (WM) (b) Figure 2: (a) Comparison of exact and approximate solutions for problem 2. (b) The variation of error norms vs t for problem 2. 5.3. Problem 3 In the third problem, we analyze Eq.(1), with parameters λ1 = λ4 = λ5 = λ6 = 0, λ2 = λ3 = λ7 = 1, σ = 1 2 , the kernel function K(t, x) = t − x, ϕ(t) = et, and the source term f(t) = te− 1 2 + e− 1 2 , and exact solution u(t) = et. The errors, Lin and Lrms for both schemes are summarized in Table 3, illustrating the convergence behavior of both schemes. While the WM shows high accuracy compared to GHQ, the latter is computationally more efficient. The exact and numerical solutions for this problem are presented in Fig. 3a, while Fig.3b provides a detailed comparison of Lin and Lrms for both numerical schemes. These results show that WM provides better accuracy than the GHQ. This analysis highlights the capability and robustness of both numerical schemes in solving this class of DIDEs. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 18 of 22 Table 3: The Lin and Lrms corresponding to Problem 3 using the GHM and WM. 1.3 GHM WM t Lin Lrms Lin Lrms 0.1 4.23×10−05 9.46×10−06 0 0 0.2 3.68×10−11 8.23×10−12 2.22×10−16 2.22×10−17 0.3 6.76×10−12 1.51×10−12 0 0 0.4 2.74×10−11 6.13×10−12 6.66×10−16 6.66×10−17 0.5 8.85×10−11 1.98×10−11 2.22×10−16 2.22×10−17 0.6 2.59×10−10 5.80×10−11 2.22×10−16 2.22×10−17 0.7 7.01×10−10 1.57×10−10 4.44×10−16 4.44×10−17 0.8 1.78×10−09 3.98×10−10 0 0 0.9 4.28×10−09 9.57×10−10 1.33×10−15 1.33×10−16 1 9.86×10−09 2.20×10−09 1.78×10−15 1.78×10−16 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 Exactsol numericsol (a) 0 0.2 0.4 0.6 0.8 1 10-18 10-16 10-14 10-12 10-10 10-8 10-6 10-4 L in (GHM) L in (WM) L rms (GHM) L rms (WM) (b) Figure 3: (a) Comparison of numerical and exact solutions for problem 3. (b) The variation of error norms vs t for problem 3. 5.4. Problem 4 In the forth problem, we examine Eq.(1), with all parameters set to λ1 = λ2 = λ3 = λ4 = λ5 = λ6 = λ7 = 1, σ = 1, the kernel function defined as K(t, x) = t − x, ϕ(t) = et, and the source term given by f(t) = t(1 + e−1) − et − et−1 + 1 + e−1. The exact solution of this problem is u(t) = et. The problem is solved using the two numerical schemes, the GHQ and WM. Table 4 provides the error measures Lin and Lrms for both schemes, which confirm their convergence. As observed in the previous problems, the WM achieves high accuracy compared to the GHQ, while the GHQ demonstrates better computational efficiency. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 19 of 22 The exact and numerical solutions for this problem are depicted in Fig. 4a, and a com- parative analysis of Lin and Lrms for the two numerical schemes is displayed in Fig. 4b. These results further illustrate the superior accuracy of the WM. This problem exemplifies the effectiveness and versatility of the proposed numerical methods for such DIDEs. Table 4: The Lin and Lrms corresponding to Problem 4 using the GHM and WM. 1.3 GHM WM t Lin Lrms Lin Lrms 0.1 1.98×10−13 4.43×10−14 0 0 0.2 1.50×10−12 3.37×10−13 2.22×10−16 2.22×10−17 0.3 7.14×10−12 1.63×10−12 2.22×10−16 3.14×10−17 0.4 2.71×10−11 6.28×10−12 2.22×10−16 3.14×10−17 0.5 8.85×10−11 2.08×10−11 2.22×10−16 3.85×10−17 0.6 2.59×10−10 6.16×10−11 2.22×10−16 3.85×10−17 0.7 7.01×10−10 1.68×10−10 2.22×10−16 3.85×10−17 0.8 1.78×10−09 4.32×10−10 2.22×10−16 3.85×10−17 0.9 4.28×10−09 1.05×10−09 4.44×10−16 5.87×10−17 1 9.86×10−09 2.44×10−09 4.44×10−16 7.36×10−17 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 Exactsol numericsol (a) 0 0.2 0.4 0.6 0.8 1 10-18 10-16 10-14 10-12 10-10 10-8 L in (GHM) L in (WM) L rms (GHM) L rms (WM) (b) Figure 4: (a) Comparison of exact and numerical solutions for problem 4. (b) The variation of error norms vs t for problem 4. 6. Conclusion In this article, a numerical technique based on the LT was developed for the numer- ical solution of DIDEs. The suggested technique first transforms the given DIDEs into algebraic equations in the Laplace domain, which are then solved for the unknown. Two Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 20 of 22 inversion techniques, the Gauss-Hermite quadrature method and the Weeks method, are then used to transform the solution back into the time domain. The efficiency and ac- curacy of the approach were validated through various examples, with comparisons of the results obtained using the two inversion methods provided via tables and figures. The Weeks method demonstrated superior performance in terms of accuracy; however, its com- putational time is higher compared to the Gauss-Hermite quadrature method. Additionally, our results demonstrated remarkable accuracy improvements when compared to those published in the literature. These findings confirm the effectiveness of the methods and provide insights into the trade-offs between the two inversion techniques. Credit authorship contribution statement K: Methodology, Writing original draft, Software, Supervision, N.J. validated the results, Conceptualization M.I.K: Writing original draft, Investigation, Conceptualiza- tion. A.A. validated the results, Conceptualization, Software N.M: Supervision, Review manuscript, Editing F.H: validated the results, Review manuscript, Editing. Acknowledgements The authors A. Aloqaily, N. Mlaiki and F. Hasan would like to thank Prince Sultan University for paying the APC and for the support through the TAS research lab. Availability of Data and Materials All the data produced or examined in this study are provided within this article. Declarations The authors state that there are no conflicts of interest. References [1] Jim M Cushing. Integrodifferential equations and delay models in population dynam- ics, volume 20. Springer Science & Business Media, 2013. [2] Cemil Tunç and Osman Tunç. On behaviours of functional volterra integro- differential equations with multiple time lags. Journal of Taibah University for Sci- ence, 12(2):173–179, 2018. [3] Erkan Cimen and Sabahattin Yatar. Numerical solution of volterra integro-differential equation with delay. J. Math. Comput. Sci, 20(3):255–263, 2020. [4] Arsalang Tang. Analysis and numerics of delay Volterra integro-differential equations. The University of Manchester (United Kingdom), 1996. [5] Zdzis law Jackiewicz. The numerical solution of volterra functional differential equa- tions of neutral type. SIAM Journal on Numerical Analysis, 18(4):615–626, 1981. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 21 of 22 [6] L Mansouri and Z Azimzadeh. Numerical solution of fractional delay volterra integro- differential equations by bernstein polynomials. Mathematical sciences, 17(4):455– 466, 2023. [7] WH Enright and Min Hu. Continuous runge-kutta methods for neutral volterra integro-differential equations with delay. Applied Numerical Mathematics, 24(2- 3):175–190, 1997. [8] LI Zaidan. Solving linear delay volterra integro-differential equations by using galerkin’s method with bernstien polynomial. J. Bablyon Appl. Sci, 20(5):1305–1313, 2012. [9] Wanyuan Ming and Chengming Huang. Collocation methods for volterra functional integral equations with non-vanishing delays. Applied mathematics and computation, 296:198–214, 2017. [10] Hermann Brunner. Recent advances in the numerical analysis of volterra functional differential equations with variable delays. Journal of computational and applied mathematics, 228(2):524–537, 2009. [11] Azzeddine Bellour and Mahmoud Bousselsal. Numerical solution of delay integro- differential equations by using taylor collocation method. Mathematical Methods in the Applied Sciences, 37(10):1491–1506, 2014. [12] Shifeng Wu and Siqing Gan. Errors of linear multistep methods for singularly per- turbed volterra delay-integro-differential equations. Mathematics and Computers in Simulation, 79(10):3148–3159, 2009. [13] Ilhame Amirali and Hülya Acar. Stability inequalities and numerical solution for neutral volterra delay integro-differential equation. Journal of Computational and Applied Mathematics, 436:115343, 2024. [14] Fathalla A Rihan, Eid H Doha, MI Hassan, and NM2604744 Kamel. Numerical treatments for volterra delay integro-differential equations. Computational Methods in Applied Mathematics, 9(3):292–318, 2009. [15] Igor 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. [16] Fazal Haq, Kamal Shah, Ghaus ur Rahman, and Muhammad Shahzad. Numerical so- lution of fractional order smoking model via laplace adomian decomposition method. Alexandria Engineering Journal, 57(2):1061–1069, 2018. [17] Fazal Haq, Kamal Shah, Asaf Khan, Muhammad Shahzad, and Ghaus ur Rahman. Numerical solution of fractional order epidemic model of a vector born disease by laplace adomian decomposition method. Punjab University Journal of Mathematics, 49(2):12–21, 2017. [18] Muhammad Asif, Kamal Shah, Bahaaeldin Abdalla, Thabet Abdeljawad, et al. Nu- merical solution of bagley–torvik equation including atangana–baleanu derivative aris- ing in fluid mechanics. Results in Physics, 49:106468, 2023. [19] Brian Davies and Brian Martin. Numerical inversion of the laplace transform: a survey and comparison of methods. Journal of computational physics, 33(1):1–32, 1979. Kamran et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6355 22 of 22 [20] Kamran, Ujala Gul, Fahad M Alotaibi, Kamal Shah, and Thabet Abdeljawad. Com- putational approach for differential equations with local and nonlocal fractional-order differential operators. Journal of Mathematics, 2023(1):6542787, 2023. [21] JAC Weideman. Gauss–hermite quadrature for the bromwich integral. SIAM Journal on Numerical Analysis, 57(5):2200–2216, 2019. [22] Jacob Andre C Weideman. Algorithms for parameter selection in the weeks method for inverting the laplace transform. SIAM Journal on Scientific Comput- ing, 21(1):111–128, 1999. [23] Kamal Shah, Muhammad Sher, Muhammad Sarwar, and Thabet Abdeljawad. Anal- ysis of a nonlinear problem involving discrete and proportional delay with application to houseflies model. AIMS Mathematics, 9(3):7321–7339, 2024. [24] W Barrett. Convergence properties of gaussian quadrature formulae. The Computer Journal, 3(4):272–277, 1961. [25] Hidetosi Takahasi and Masatake Mori. Estimation of errors in the numerical quadra- ture of analytic functions. Applicable Analysis, 1(3):201–229, 1971. [26] Arthur H Stroud, Don Secrest, et al. Gaussian quadrature formulas, volume 374. Prentice-Hall Englewood Cliffs, NJ, 1966. [27] Joseph Abate, Gagan L Choudhury, and Ward Whitt. Numerical inversion of mul- tidimensional laplace transforms by the laguerre method. Performance Evaluation, 31(3-4):229–243, 1998. [28] Nur Inshirah Naqiah Ismail and Zanariah Abdul Majid. Numerical solution on neutral delay volterra integro-differential equation. Bulletin of the Malaysian Mathematical Sciences Society, 47(3):85, 2024. [29] Nur Inshirah Naqiah Ismail, Zanariah Abdul Majid, Norazak Senu, and Nadihah Wahi. Numerical method in solving neutral and retarded volterra delay integro- differential equations. In Journal of Physics: Conference Series, volume 1988, page 012033. IOP Publishing, 2021. [30] Aman Jhinga, Jayvant Patade, and Varsha Daftardar-Gejji. Solving volterra integro- differential equations involving delay: a new higher order numerical method. arXiv preprint arXiv:2009.11571, 2020.