EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 5766 ISSN 1307-5543 – ejpam.com Published by New York Business Global Convergence Analysis of Multi-Step Collocation Method to First-Order Volterra Integro-Differential Equation with Non-Vanishing Delay Ahmed Ali Eashel1, Saeed Pishbin1,∗, Parviz Darania1 1 Department of Mathematics, Faculty of Science, Urmia University, P.O.Box 165, Urmia- Iran Abstract. Generally, solutions to functional equations involving non-vanishing delays tend to exhibit lower regularity compared to those of smooth functions. In this context, we examine a first- order Volterra integro-differential equation (VIDE) with a non-vanishing delay, delving into the characteristics of its solutions. To enhance the accuracy of traditional one-step collocation methods [1], we employ multi-step collocation techniques to obtain numerical solutions for the VIDE with non-vanishing delay. The global convergence properties of the multi-step numerical approach are scrutinized using the Peano Kernel Theorem. Subsequently, for comparative analysis, we utilize a one-step collocation method to numerically solve this equation, showcasing the effectiveness and precision of the multi-step collocation method. 2020 Mathematics Subject Classifications: 65R20, 65L03 Key Words and Phrases: Volterra integro-differential equation, Delay integro-differential equa- tion, Multi-step collocation methods, Convergence analysis 1. Introduction Here, we deal with the approximate solution of the delay first-order Volterra integro- differential equation of the form x′(t) = c1(t)x(t)+c2(t)x(τ(t))+f(t)+ ∫ t t0 K(t, s, x(s))ds+ ∫ τ(t) t0 K̂(t, s, x(s))ds, t ∈ J = [t0, T ], (1) where x(t) = ζ(t), t ∈ [τ(t0), t0], x(t) is the unknown solution and c1, c2,K, K̂ are given functions. Also τ(t) is delay (or lag) function which will be defined completely in section 2. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5766 Email addresses: a.alieashel@gmail.com (A. Ali Eashel), s.pishbin@urmia.ac.ir (S. Pishbin), p.darania@urmia.ac.ir (P. Darania) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 2 of 29 Volterra integro-differential equations (VIDEs) may be viewed as ordinary differential equations affected by a memory term arising from an integral operator which exhibit a common characteristic, when the initial data was smooth, it resulted in smooth solutions. In general this statement does not hold for an equation involving a non-vanishing delay. In such cases, the presence of delays gives rise to primary points of discontinuity, resulting in a solution that is less regular compared to the initially supplied smooth functions. The initial-value problem for Volterra integro-differential equation arise in some math- ematical modelling processes in biological and physical phenomena, such as fluid dynamics, viscoelasticity in materials with memory, such as population dynamics [2–6]. Also, VIDEs with non-vanishing delay have been used in describing some phenomena of population growth the transmission of an epidemic with the influx of immigrants into the population, incorporating scenarios into models like the predator-prey model. For instance, in [7] , we can find the following delay integro-differential equations (IDEs) related to population dynamics: z′(t) = cz(t) ( 1− 1 L ∫ 0 −τ z(t+ s)dσ(s) ) , where L is the environmental carrying capacity and c is the intrinsic rate. Also, VIDEs with non-vanishing delay have applications in economics and Belairs model to life span (See examples on pages 212-214 in [1] for more details). For the numerical soltion of Volterra integro-differential equations, we can find several numerical methods in the literature, i.e. rationalized Haar functions [8], Galerkin and collocation type methods [1, 9–12], Runge-Kutta methods [1, 9], Lagrange polynomial method [13] and multi-step methods [14]. In [14], the authors, have analyzed multi-step collocation methods for Volterra IDEs and derived order of convergence of the proposed method. The study involved an examination of numerical stability analysis, and various classes of methods that are A0-stable were presented. In [15], the authors considered the Chebyshev collocation method to solve the IDEs and developed this method for the system of Volterra-Fredholm IDEs. Also, polynomial approximation based on Taylor expansion [16] has been developed for the Volterra-Fredholm and VIDEs in real application. Integro- differential equations involving convolution integrals as a generalization of the fractional differential equations has been solved by numerical method in [17] and has been shown that the proposed numerical scheme is stable. The authors, [18] consider a graded mesh refinement algorithm for solving time-delayed parabolic partial differential equations with a small diffusion parameter. Also, in [19] a class of boundary layer originated singularly perturbed parabolic reaction-diffusion problems with large time delay have been studied. A nonlinear system of singularly perturbed delay differential equation whose each com- ponent of the solution has multiple layers has been investigated in [20]. Delay integral equations (IEs) have been solved approximately by many authors (see, e.g.,[1, 21–28]). Recently, some authors are interested in working on the numerical solution of delay IDEs. In [29], the authors used multi-step methods to find a approximate solution of singularly perturbed delay VIDEs. One- step polynomial collocation method [1] has been applied to find numerical solution of delay VIDEs by Brunner. He performed the convergence analysis of the collocation method using Peano kernel theorem and investigated super convergence A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 3 of 29 analysis of the numerical approximation in detail. The delay VIDEs have been achieved through the utilization of the two-point multi-step block (2PBM) [30] method which was formulated by Taylor expansion. The authors developed the 2PBM method by consider- ing the predictor-corrector formulae. Also, in [31] an effective numerical method has been applied to solve VIDEs including neutral terms with variable delays by fundamental ma- trices of Laguerre polynomials. Here, we consider multi-step collocation methods to delay equation (1) and try to increase the order of convergence in comparing of the one-step collocation methods in [1]. The paper is organized as follows: Section 2 is dedicated to proposing the structure of solutions for VIDEs with non- vanishing delays and subsequently outlining the application of the multi-step polynomial collocation method by employing well-known interpolation polynomials. In Section 3, we construct and analyze the convergence of numerical solutions, taking into consideration Peano’s Theorem for interpolation. This is succeeded by the discussion of two test prob- lems in Section 4 to validate the theoretical results. Finally, in Section 5, we conclude the paper and suggest potential future avenues for research, which are currently less explored. 2. Numerical method for the delay VIDE In this section, firstly we will state structure of the solution of VIDEs with non- vanishing delays (1) and then consider the multi-step collocation method to solve this equation numerically. 2.1. Structure of the solution of VIDEs with non-vanishing delays Let τ(t) = t − α(t) be strictly increasing on J with α(t) ≥ α0 > 0 for t ∈ J and α(t) ∈ Cν(J) for some ν ≥ 0. The presence of a non-vanishing delay τ(t) gives rise to the primary discontinuity points denoted as ςi so that they are obtained from the following formula τ(ςi) = ςi−1, i ≥ 1, ς0 = t0, and ςµ − ςµ−1 = α(ςµ) ≥ α0 > 0, for all µ ≥ 0. Now, we write the equation (1) in the local form x′(t) = c1(t)x(t) + q1,i(t) + ∫ t ςi K(t, s, x(s))ds, t ∈ [ςi, ςi+1], (2) where q1,i(t) = c2(t)x(τ(t)) + f(t) + ∫ ςi t0 K(t, s, x(s))ds+ ∫ τ(t) t0 K̂(t, s, x(s))ds, i ≥ 1. For t ∈ (ς0, ς1], we derive q1,0(t) = f(t) + c2(t)x(τ(t))− ∫ t0 τ(t) K̂(t, s, ζ(s))ds. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 4 of 29 and lim t→t−0 x′(t) = ζ ′(t0), lim t→t+0 x′(t) = c1(t0)x(t0) + q1,0(t0), then, x′ has a discontinuity at t = t0. Assume that c1, c2, f ∈ C(J), K(., .) ∈ C(D × R) and K̂(., .) ∈ C(Dτ × R), with D = {(t, s) : t0 ≤ s ≤ t ≤ T}, Dτ = {(t, s) : τ(t0) ≤ s ≤ τ(t), t ∈ J}. For t ∈ (ς1, ς2], we have x′(ς−1 ) = c1(ς1)x(ς1) + q1,0(ς − 1 ) + ∫ ς1 ς0 K(ς1, s, x(s))ds, x′(ς+1 ) = c1(ς1)x(ς1) + q1,1(ς + 1 ). Then x′(ς−1 )− x′(ς+1 ) = 0, and for t ∈ [ςi, ςi+1], i ≥ 2, this continuity is maintained for the other points of ςi as well. Now, using these arguments and a similar process in the proof of Theorem 2.2 from [32], we consider the following theorem which gives the relevant conditions for the investigation of the unique solution of the equation (1): Theorem 1. Assume that τ(t) = t− α(t) be strictly increasing on J with α(t) ≥ α0 > 0 for t ∈ J and α(t) ∈ Cν(J) for some ν ≥ 0. Also 1. c1, c2, f ∈ C(J) and C1 = max t∈J |c1(t)|. 2. K(., .) ∈ C(D × R) and K̂(., .) ∈ C(Dτ × R) with D = {(t, s) : t0 ≤ s ≤ t ≤ T}, Dτ = {(t, s) : τ(t0) ≤ s ≤ τ(t), t ∈ J}. 3. K satisfies the Lipschitz condition |K(t, s, x)−K(t, s, y)| ≤ L|x− y| ∀(t, s) ∈ D,x, y ∈ R. So it can be said for each initial function ζ(t) ∈ C[τ(t0), t0] the equation (1) has a unique solution x ∈ C(J) ∩ C1(t0, T ]. Also, in general, at t = t0 its derivative is discontinuous (but bounded): lim t→t−0 x′(t) ̸= lim t→t+0 x′(t). Proof. We should consider local form of the equation (1) by (2). For µ = 0, we have t ∈ [ς0, ς1] and derive x′(t) = c1(t)x(t) + q1,0(t) + ∫ t t0 K(t, s, x(s))ds, x(t0) = ζ(t0). (3) A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 5 of 29 It is equivalent to a nonlinear VIE of the second kind as: x(t) = q̂1,0(t) + ∫ t t0 c1(s)x(s)ds+ ∫ t t0 ∫ t s K(η, s, x(s))dηds, (4) where q̂1,0(t) = ζ(t0) + ∫ t t0 q1,0(s)ds. We rewrite equation (4) in operator form x(t) = q̂1,0(t) + V(x)(t), (5) where V(x)(t) = ∫ t t0 c1(s)x(s)ds+ ∫ t t0 ∫ t s K(η, s, x(s))dηds. We choose a positive constant δ0 > 0 and define A1 = C([t0, t0 + δ0]). Then A1 is a Banach space equipped with the maximum norm. We will prove that the operator V restricted to A1 is a contraction operator. For x, y ∈ A1 |V(x)− V(y)| ≤ ∫ t t0 |c1(s)| |x(s)− y(s)|ds+ ∫ t t0 ∫ t s |K(η, s, x(s))−K(η, s, y(s))|dηds ≤ C1δ0∥x− y∥∞ + Lδ20 2 ∥x− y∥∞ ≤ ∥x− y∥∞β(δ0 + δ20), (6) where β = max{C1, L 2 }. We know that β > 0, then, for 0 < δ0 < −1+ √ 1+ 4 β 2 , V is a contraction map on the Banach space A1, and has a unique fixed point x1 ∈ A1. In the sequel, we study the case t ∈ [t0 + δ0, t0 + 2δ0] if t0 + 2δ0 < ς1. By defining A2 = {x ∈ C([t0, t0 + 2δ0]), and x(t) = x1(t), for t ∈ [t0, t0 + δ0]}, we will show that V is also a contraction map on A2: |V(x)− V(y)| ≤ ∫ t t0 |c1(s)| |x(s)− y(s)|ds+ ∫ t t0 ∫ t s |K(η, s, x(s))−K(η, s, y(s))|dηds ≤ ∫ t t0 |c1(s)| |x(s)− y(s)|ds+ ∫ t t0 (t− s)L|x(s)− y(s)|ds ≤ ∫ t0+δ0 t0 |c1(s)| |x1(s)− x1(s)|ds+ ∫ t0+δ0 t0 (t− s)L|x1(s)− x1(s)|ds + ∫ t t0+δ0 |c1(s)| |x(s)− y(s)|ds+ ∫ t t0+δ0 (t− s)L|x(s)− y(s)|ds ≤ C1δ0∥x− y∥∞ + Lδ20 2 ∥x− y∥∞ ≤ ∥x− y∥∞β(δ0 + δ20). (7) For 0 < δ0 < −1+ √ 1+ 4 β 2 , it also follows from the Banach fixed point theorem that V has a unique fixed point x2 ∈ A2, which satisfies x2(t) = x1(t) for t ∈ [t0, t0 + δ0]. Therefore, we obtain a unique solution to VIDE (3) on the interval [t0, t0 + 2δ0]. In the sequel, we A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 6 of 29 proceed similarly to the strategy considered in the proof of the Theorem 2.2 from [32] and complete the proof. Uniqueness: Suppose that (4) possesses two continuous solutions x and z on the interval [ς0, ς1] . Hence, by the Lipschitz condition: |x(t)− z(t)| ≤ ∫ t t0 |c1(s)| |x(s)− z(s)|ds+ ∫ t t0 ∫ t s |K(η, s, x(s))−K(η, s, z(s))|dηds ≤ (C1 + L(ς1 − ς0)) ∫ t t0 |x(s)− z(s)|ds. (8) It follows from the classical Gronwall lemma [1] |x(t)− z(t)| ≤ 0× Exp[(C1 + L(ς1 − ς0))(t− t0)] = 0, t ∈ [ς0, ς1]. (9) The continuity of x and z then implies that x(t) = z(t) for all t ∈ [ς0, ς1]. For µ ≥ 1, the above discussion (µ = 0) is readily adapted to establish the (local) existence and uniqueness of a solution. If the data in the delay VIDE (1) are smooth functions, the corresponding solution will essentially inherit this smoothness, except at the primary discontinuity points ςµ. Considering similar strategy considered in the Theorem 2.3 from [32] for the local form (2), we derive a regularity result for the solution of the equation (1) as: Theorem 2. Assume that τ(t) = t− α(t) be strictly increasing on J with α(t) ≥ α0 > 0 for t ∈ J and α(t) ∈ Cν(J) for some ν ≥ d. Also 1. c1, c2, f ∈ Cd(J) and ζ(t) ∈ Cd[τ(t0), t0]. 2. K(., .) ∈ Cd(D ×R) and K̂(., .) ∈ Cd(Dτ ×R). 3. K satisfies the Lipschitz condition |K(t, s, x)−K(t, s, y)| ≤ L|x− y| ∀(t, s) ∈ D,x, y ∈ R. The unique solution of the equation (1) is (d+ 1)-times continuously differentiable on each left-open macro-interval (ςµ, ςµ+1] for each µ = 0, 1, . . . ,M and has a bounded first derivative on J . Also, in general, at t = ςµ, (µ = 0, 1, . . . ,min{d,M}), we have: lim t→ς−µ x(µ)(t) = lim t→ς+µ x(µ)(t), while the (µ+ 1)st derivative of x is in general not continuous at t = ςµ. If min{d,M} = d < M , the solution possesses a continuous (d+ 1)st derivative on [ςµ, T ]. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 7 of 29 Remark 1. If the data in the delay VIDE (1) are smooth functions with the degree of smoothness d, the corresponding solution will essentially inherit this smoothness on each left-open macro-interval (ςµ, ςµ+1]. Also, for t = ςµ(µ = 0, 1, · · · ,min{d,M}), the µ st derivative of the solution is continuous at this points. Since solutions of this problem generally suffer from a loss of regularity at the primary discontinuity points ςµ, the mesh lh = {tn : t0 < t1 < · · · < tN} underlying the collocation space will have to include these points if the collocation solution is to attain its optimal global or local order of convergence. Thus, we shall employ meshes of the form Lh := M⋃ µ=0 l (µ) h , l (µ) h := {t(µ)n : ςµ = t (µ) 0 < t (µ) 1 < · · · < t (µ) Nµ = ςµ+1}. Such a mesh is called a constrained mesh (with respect to τ(t)) for J . We will refer to Lh as the macro-mesh and call the l (µ) h the underlying local meshes. Maybe we can consider the other ideas [33–35] to overcome the problem of low smoothness of the solution. 2.2. Multi-step method In this subsection, we apply the multi-step collocation method to solve the equation (1). Let T in J = [t0, T ] is defined so that T = ςM+1, for some M ≥ 1, Lh := M⋃ µ=0 l (µ) h , l (µ) h := {t(µ)n : ςµ = t (µ) 0 < t (µ) 1 < · · · < t (µ) Nµ = ςµ+1}, h (µ) n = t (µ) n+1 − t (µ) n , µ = 0, . . . ,M, (M ≥ 1) and for 0 ≤ n ≤ N − 1, Y (µ) h = {t(µ)n,i = t(µ)n + sih (µ) n : 0 < s1 < · · · < sm ≤ 1}, where {si} are collocation parameters. We consider w′ as an approximate solution of x′ in [t (µ) n , t (µ) n+1] by w′(t(µ)n + zh(µ)n ) = r−1∑ k=0 Pk(z)x ′(µ) n−k + m∑ j=1 Qj(z)W (µ) n,j , z ∈ (0, 1], W (µ) n,j = w′(t (µ) n,j), n ≥ r − 1, (10) where x ′(µ) n−k = w′(t (µ) n−k) and Pk(z) = ( m∏ i=1 z − si −k − si )( r−1∏ i=0,i ̸=k z + i −k + i ) , Qj(z) = ( r−1∏ i=0 z + i sj + i )( m∏ i=1,i ̸=j z − si sj − si ) . (11) A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 8 of 29 Now, setting x (µ) n = w(t (µ) n ) and αk(z) = ∫ z 0 Pk(s)ds, (k = 0, · · · , r − 1), βj(z) = ∫ z 0 Qj(s)ds, (j = 1, · · · ,m), we obtain from (10) w(t(µ)n + zh(µ)n ) = x(µ)n + h(µ)n r−1∑ k=0 αk(z)x ′(µ) n−k + h(µ)n m∑ j=1 βj(z)W (µ) n,j , (12) and x ′(µ) n+1 = r−1∑ k=0 Pk(1)x ′(µ) n−k + m∑ j=1 Qj(1)W (µ) n,j , x (µ) n+1 = x(µ)n + h(µ)n r−1∑ k=0 αk(1)x ′(µ) n−k + h(µ)n m∑ j=1 βj(1)W (µ) n,j , n ≥ r − 1. In [t (µ) n , t (µ) n+1], 0 ≤ n < r−1, the primary values x ′(µ) 0 , x ′(µ) 1 , x ′(µ) 2 , . . . , x ′(µ) r−1 and x (µ) r−1 may be obtained using the appropriate methods (see one- step collocation method in chapter 3 in [1]). Also, when t = t0, we have x′(t0) = c1(t0)x(t0) + q1,0(t0) where x(t0) = ζ(t0). In addition, w as an approximate solution should satisfy the following collocation equation w′(t (µ) n,i ) = c1(t (µ) n,i )w(t (µ) n,i ) + c2(t (µ) n,i )w(τ(t (µ) n,i )) + f(t (µ) n,i ) + ∫ t (µ) n,i t0 K(t (µ) n,i , s, w(s))ds + ∫ τ(t (µ) n,i ) t0 K̂(t (µ) n,i , s, w(s))ds. (13) Let τ(t) be linear. Inserting (10),(12) into (13) and using appropriate change of vari- ables for each sub interval [t (µ) n , t (µ) n+1], we have the following non-linear system in two cases: I) For µ = 0, we have W (0) n,i = c1(t (0) n,i) ( x (0) n + h (0) n r−1∑ k=0 αk(si)x ′(0) n−k + h(0)n m∑ j=1 βj(si)W (0) n,j ) + f(t (0) n,i) + r−2∑ l=0 h (0) l ∫ 1 0 K ( t (0) n,i , t (0) l + sh (0) l , w(t (0) l + sh (0) l ) ) ds + n−1∑ l=r−1 h (0) l ∫ 1 0 K ( t (0) n,i , t (0) l + sh (0) l , x (0) l + h (0) l r−1∑ k=0 αk(s)x ′(0) l−k + h (0) l m∑ j=1 βj(s)W (0) l,j ) ds +h (0) n ∫ si 0 K ( t (0) n,i , t (0) n + sh(0)n ), x(0)n + h(0)n r−1∑ k=0 αk(s)x ′(0) n−k + h(0)n m∑ j=1 βj(s)W (0) n,j ) ds +c2(t (0) n,i)ζ(τ(t (0) n,i)) + ∫ τ(t (0) n,i) t0 K̂(t (0) n,i , s, ζ(s))ds. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 9 of 29 II) For µ ≥ 1 W (µ) n,i = c1(t (µ) n,i ) ( x (µ) n + h (µ) n r−1∑ k=0 αk(si)x ′(µ) n−k + h(µ)n m∑ j=1 βj(si)W (µ) n,j ) + f(t (µ) n,i ) + µ−1∑ η=0 r−2∑ l=0 h (η) l ∫ 1 0 K ( t (µ) n,i , t (η) l + sh (η) l , w(t (η) l + sh (η) l ) ) ds + µ−1∑ η=0 Nµ−1∑ l=r−1 h (η) l ∫ 1 0 K ( t (µ) n,i , t (η) l + sh (η) l , x (η) l + h (η) l r−1∑ k=0 αk(s)x ′(η) l−k + h (η) l m∑ j=1 βj(s)W (η) l,j ) ds + r−2∑ l=0 h (µ) l ∫ 1 0 K ( t (µ) n,i , t (µ) l + sh (µ) l , w(t (µ) l + sh (µ) l ) ) ds + n−1∑ l=r−1 h (µ) l ∫ 1 0 K ( t (µ) n,i , t (µ) l + sh (µ) l , x (µ) l + h (µ) l r−1∑ k=0 αk(s)x ′(µ) l−k + h (µ) l m∑ j=1 βj(s)W (µ) l,j ) ds +h (µ) n ∫ si 0 K ( t (µ) n,i , t (µ) n + sh(µ)n ), x(µ)n + h(µ)n r−1∑ k=0 αk(s)x ′(µ) n−k + h(µ)n m∑ j=1 βj(s)W (µ) n,j ) ds +c2(t (µ) n,i ) ( x (µ−1) n + h (µ−1) n r−1∑ k=0 αk(si)x ′(µ−1) n−k + h(µ−1) n m∑ j=1 βj(si)W (µ−1) n,j ) + µ−2∑ η=0 r−2∑ l=0 h (η) l ∫ 1 0 K̂ ( t (µ) n,i , t (η) l + sh (η) l , w(t (η) l + sh (η) l ) ) ds + µ−2∑ η=0 Nµ−1∑ l=r−1 h (η) l ∫ 1 0 K̂ ( t (µ) n,i , t (η) l + sh (η) l , x (η) l + h (η) l r−1∑ k=0 αk(s)x ′(η) l−k + h (η) l m∑ j=1 βj(s)W (η) l,j ) ds + r−2∑ l=0 h (µ−1) l ∫ 1 0 K̂ ( t (µ) n,i , t (µ−1) l + sh (µ−1) l , w(t (µ−1) l + sh (µ−1) l ) ) ds + n−1∑ l=r−1 h (µ−1) l ∫ 1 0 K̂ ( t (µ) n,i , t (µ−1) l + sh (µ−1) l , x (µ−1) l + h (µ−1) l r−1∑ k=0 αk(s)x ′(µ−1) l−k +h (µ−1) l m∑ j=1 βj(s)W (µ−1) l,j ) ds +h (µ−1) n ∫ si 0 K̂ ( t (µ) n,i , t (µ−1) n + sh(µ−1) n ), x(µ−1) n + h(µ−1) n r−1∑ k=0 αk(s)x ′(µ−1) n−k +h (µ−1) n m∑ j=1 βj(s)W (µ−1) n,j ) ds. We can get WWW (µ) n by solving the non-linear system. Substituting it into (12), the approximate solution of (1) can be obtained. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 10 of 29 Remark 2. If you consider linear case of the equation (1) as: x′(t) = c1(t)x(t)+c2(t)x(τ(t))+f(t)+ ∫ t t0 K(t, s)x(s)ds+ ∫ τ(t) t0 K̂(t, s)x(s)ds, t ∈ J = [t0, T ], (14) where x(t) = ζ(t), t ∈ [τ(t0), t0]. Then we have the following linear algebraic system for µ = 0 and µ ≥ 1 as: I) For µ = 0, we have W (0) n,i = c1(t (0) n,i) ( x (0) n + h (0) n r−1∑ k=0 αk(si)x ′(0) n−k + h(0)n m∑ j=1 βj(si)W (0) n,j ) + f(t (0) n,i) + r−2∑ l=0 h (0) l ∫ 1 0 K(t (0) n,i , t (0) l + sh (0) l )w(t (0) l + sh (0) l )ds + n−1∑ l=r−1 h (0) l ∫ 1 0 K(t (0) n,i , t (0) l + sh (0) l ) ( x (0) l + h (0) l r−1∑ k=0 αk(s)x ′(0) l−k + h (0) l m∑ j=1 βj(s)W (0) l,j ) ds +h (0) n ∫ si 0 K(t (0) n,i , t (0) n + sh(0)n ) ( x(0)n + h(0)n r−1∑ k=0 αk(s)x ′(0) n−k + h(0)n m∑ j=1 βj(s)W (0) n,j ) ds +c2(t (0) n,i)ζ(τ(t (0) n,i)) + ∫ τ(t (0) n,i) t0 K̂(t (0) n,i , s)ζ(s)ds. II) For µ ≥ 1 A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 11 of 29 ( Im − h (µ) n (CCC (µ) 1n βββ + h (µ) n DDDn,µ n ) ) WWW (µ) n = x (µ) n CCC (µ) 1n + h (µ) n CCC (µ) 1n αααXXX (µ) n +FFF (µ) n + µ−1∑ η=0 r−2∑ l=0 h (η) l ZZZ l,η n + µ−1∑ η=0 Nµ−1∑ l=r−1 h (η) l (x (η) l VVV l,η n +h (η) l GGGl,η n XXX (η) l + h (η) l DDDl,η n WWW (η) l ) + r−2∑ l=0 h (µ) l ZZZ l,µ n + n−1∑ l=r−1 h (µ) l (x (µ) l VVV l,µ n + h (µ) l GGGl,µ n XXX (µ) l +h (µ) l DDDl,µ n WWW (µ) l ) + h (µ) n (x (µ) n VVV n,µ n + h (µ) n GGGn,µ n XXX (µ) n ) +x (µ−1) n CCC (µ−1) 2n + h (µ−1) n CCC (µ−1) 2n αααXXX (µ−1) n + h (µ−1) n CCC (µ−1) 2n βββ + µ−2∑ η=0 r−2∑ l=0 h (η) l ẐZZ l,η n + r−2∑ l=0 h (µ−1) l ẐZZ l,µ−1 n + µ−2∑ η=0 Nµ−1∑ l=r−1 h (η) l (x (η) l V̂VV l,η n + h (η) l ĜGG l,η n XXX (η) l + h (η) l D̂DD l,η n WWW (η) l ) + n−1∑ l=r−1 h (µ−1) l (x (µ−1) l V̂VV l,µ−1 n + h (µ−1) l ĜGG l,µ−1 n XXX (µ−1) l +h (µ−1) l D̂DD l,µ−1 n WWW (µ−1) l ) + h (µ−1) n (x (µ−1) n V̂VV n,µ−1 n +h (µ−1) n ĜGG n,µ−1 n XXX (µ−1) n + h (µ−1) n D̂DD n,µ−1 n ), (15) where βββ = ( βj(si) i, j = 1, · · · ,m ) , ααα =  αk(si) i,= 1, . . . ,m, k = 0, . . . , r − 1,  , WWW (.) l = ( W (.) l,1 , . . . ,W (.) l,m )T , CCC (.) 1n = diag ( c1(t (.) n,1), . . . , c1(t (.) n,m) ) , CCC (.) 2n = diag ( c2(t (.) n,1), . . . , c2(t (.) n,m) ) , ( DDDl,. n ) i,j =  ∫ 1 0 K(t (µ) n,i , t (.) l + sh (.) l )βj(s)ds, l = r − 1, . . . , n− 1, ∫ si 0 K(t (µ) n,i , t (.) l + sh (.) l )βj(s)ds, l = n, i, j = 1, . . . ,m, A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 12 of 29 ( GGGl,. n ) i,k =  ∫ 1 0 K(t (µ) n,i , t (.) l + sh (.) l )αk(s)ds, l = r − 1, · · · , n− 1, ∫ si 0 K(t (µ) n,i , t (.) l + sh (.) l )αk(s)ds, l = n, i = 1, . . . ,m, k = 0, . . . , r − 1, XXX (.) l = (x ′(.) l , · · · , x′(.)l−r+1) T , FFF (µ) n = (f(t (µ) n,1), · · · , f(t (µ) n,m))T , ZZZ l,. n = (∫ 1 0 K(t (µ) n,1, t (.) l +sh (.) l )w(t (.) l +sh (.) l )ds, . . . , ∫ 1 0 K(t(µ)n,m, t (.) l +sh (.) l )w(t (.) l +sh (.) l )ds )T , VVV l,. n =  (∫ 1 0 K(t (µ) n,1, t (.) l + sh (.) l )ds, . . . , ∫ 1 0 K(t(µ)n,m, t (.) l + sh (.) l )ds )T , l = r − 1, . . . , n− 1, (∫ s1 0 K(t (µ) n,1, t (.) n + sh(.)n )ds, . . . , ∫ sm 0 K(t(µ)n,m, t(.)n + sh(.)n )ds )T , l = n. Also the matrices D̂DD l,. n , ĜGG l,. n , ẐZZ l,. n and V̂VV l,. n are defined similarly to the above matrices, only instead of K, we put K̂. By solving the linear system obtained above, we can get WWW (µ) n and substituting it into (12), the approximate solution of (14) can be achieved. 3. Convergence analysis In this section, we consider convergence analysis of the proposed numerical method for the linear case (14) and at the end of this section, we will explain how to extended it for the non-linear case. Remark 3. In [14], the authors applied multi-step collocation methods for classical VIDEs and in the Theorem 3.1, showed that the order of convergence of this method is m+ r− 1. Note that in this paper, they approximated exact solution instead of the derivative of the exact solution. If they approximate the derivative of the exact solution and by integrating get the approximate of the exact solution then they could archive the order of convergence m+ r. Theorem 3. Assume that for d ≥ m+ r, cϑ ∈ Cd(I), ϑ = 1, 2,K ∈ Cd(D) K̂ ∈ Cd(Dτ ) and ζ ∈ Cd+1([τ(t0), t0]). Let τ(t) = t−α(t), be strictly increasing on J with α(t) ≥ α0 > 0 for t ∈ J and α(t) ∈ Cd(J). Also, for PPP = [ 000r−1,1 Ir−1 Pr−1(1) Pr−2(1), · · · , P0(1), ] , A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 13 of 29 the spectral radius of the matrix, denoted by ρ(PPP ), is less than 1 and the starting errors are || ε ||∞,[t (µ) 0 ,t (µ) r−1] = O(h(µ))p. Then the estimates ∥e(γ)∥∞ = ∥x(γ) − w(γ)∥∞ ≤ Cγh p, γ = 0, 1, (16) hold for any collocation parameters {si} in [0, 1], and h = max l,ν h (ν) l , p = m+ r. Proof. Assume that e := x−w show the collocation error for the approximate solution w in (12) which satisfies in the following equation e′(t) = c1(t)e(t)+c2(t)e(τ(t))+δ(t)+ ∫ t t0 K(t, s)e(s)ds+ ∫ τ(t) t0 K̂(t, s)e(s)ds, t ∈ J, (17) e(t) = 0, t ∈ [τ(t0), t0]. Also, δ(t) = 0, t ∈ Yh = M⋃ µ=0 Y (µ) h . Considering (17), for t ∈ I(µ) = (ςµ, ςµ+1], we have e′(t) = c1(t)e(t) + Gµ(t) + δ(t) + ∫ t ςµ K(t, s)e(s)ds, t ∈ I(µ), (18) where Gµ(t) = c2(t)e(τ(t)) + ∫ ςµ t0 K(t, s)e(s)ds+ ∫ τ(t) t0 K̂(t, s)e(s)ds. (19) For µ = 0, it follows that e′(t) = c1(t)e(t) + δ(t) + ∫ t t0 K(t, s)e(s)ds, t ∈ I(0). (20) Using this remark on I(0), we can obtain error bound as: ∥ e(ν) ∥∞≤ Cν(h (0))p, ν = 0, 1. (21) and e(ν)(ς1) = O((h (0) l )p). Now, on the interval I(µ), 1 ≤ µ ≤ M , the collocation error equation (17) in t = t (µ) n,i , satisfies the equation e′(t (µ) n,i ) = c1(t (µ) n,i )e(t (µ) n,i ) + c2(t (µ) n,i )e(τ(t (µ) n,i )) + δ(t (µ) n,i ) + ∫ t (µ) n,i t0 K(t (µ) n,i , s)e(s)ds+ ∫ τ(t (µ) n,i ) t0 K̂(t (µ) n,i , s)e(s)ds, t ∈ I(µ), (22) after some computation equation (22) reduce the following form A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 14 of 29 e′(t (µ) n,i ) = δh(t (µ) n,i ) + c1(t (µ) n,i )e(t (µ) n,i ) + c2(t (µ) n,i )e(t (µ−1) n,i ) + µ−1∑ v=0 ∫ ςv+1 ςv K(t (µ) n,i , s)e(s)ds + r−1∑ l=1 h (µ) l ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )e(t (µ) l + sh (µ) l )ds + n−1∑ l=r h (µ) l ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )e(t (µ) l + sh (µ) l )ds +h (µ) n ∫ si 0 K(t (µ) n,i , t (µ) n + sh(µ)n )e(t(µ)n + sh(µ)n )ds + µ−2∑ v=0 ∫ ςv+1 ςv K̂(t (µ) n,i , s)e(s)ds + r−1∑ l=1 h (µ−1) l ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )e(t (µ−1) l + sh (µ−1) l )ds + n−1∑ l=r h (µ−1) l ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )e(t (µ−1) l + sh (µ−1) l )ds +h (µ−1) n ∫ s̃i 0 K̂(t (µ) n,i , t (µ−1) n + sh(µ−1) n )e(t(µ−1) n + sh(µ−1) n )ds. (23) By the hypothesis on the starting error it follows that ε(t (µ) l + vh (µ) l ) = (h (µ) l )m+rql(s), l = 0, . . . , r − 2, s ∈ (0, 1], (24) with ∥ql∥∞ ≤ C1 independent of h (µ) l . Recall now the analogous error equations for equation (17), for e and e′, they are, respectively, e(t(µ)n +zh(µ)n ) = e(t(µ)n )+h(µ)n r−1∑ k=0 αk(z)e ′(µ) n−k+h(µ)n m∑ j=1 βj(z)E (µ) n,j+(h(µ)n )p+1R (µ) m+r,n(z), (25) and e′(t(µ)n + zh(µ)n ) = r−1∑ k=0 Pk(z)e ′(µ) n−k + m∑ j=1 Qj(z)E (µ) n,j + (h(µ)n )pR ′(µ) m+r,n(z), z ∈ (0, 1], (26) where R ′(µ) d,n (s) = ∫ 1 1−r Kd,r(s, z)e ′(d)(t(µ)n + zh(µ)n )dz, A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 15 of 29 Kd,r(s, z) = 1 (d− 1)! (s− z)d−1 + − r−1∑ k=0 Pk(s)(−k − z)d−1 + − m∑ j=1 Qj(s)(sj − z)d−1 +  , with E (µ) n,i = X (µ) n,i −W (µ) n,i . Substituting the equations (25) and (26) in the equation (23), we get E (µ) n,i = c1(t (µ) n,i )e(t (µ) n ) + c1(t (µ) n,i )h (µ) n r−1∑ k=0 αk(si)e ′(µ) n−k + c1(t (µ) n,i )h (µ) n m∑ j=1 βj(si)E (µ) n,j +c2(t (µ) n,i )e(t (µ−1) n ) + c2(t (µ) n,i )h (µ) n r−1∑ k=0 αk(si)e ′(µ−1) n−k + c2(t (µ) n,i )h (µ) n m∑ j=1 βj(si)E (µ−1) n,j +c1(t (µ) n,i )(h (µ) n )p+1R (µ) m+r,n(si) + c2(t (µ) n,i )(h (µ−1) n )p+1R (µ−1) m+r,n(si) + Nµ−1∑ l=0 µ−1∑ v=0 (h (v) l )p+2 ∫ 1 0 K(t (µ) n,i , t (v) l + sh (v) l )R (v) m+r,l(s)ds + n−1∑ l=r (h (µ) l )p+2 ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )R (µ) m+r,l(s)ds +(h (µ) n )p+2 ∫ si 0 K(t (µ) n,i , t (µ) n + sh(µ)n )R (µ) m+r,n(s)ds + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )p+2 ∫ 1 0 K̂(t (µ) n,i , t (v) l + sh (v) l )R (v) m+r,l(s)ds + n−1∑ l=r (h (µ−1) l )p+2 ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )R (µ−1) m+r,l(s)ds +(h (µ−1) n )p+2 ∫ si 0 K̂(t (µ) n,i , t (µ−1) n + sh(µ−1) n )R (µ−1) m+r,n(s)ds + r−1∑ l=1 (h (µ) l )p+1 ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l ql(s)ds + r−1∑ l=1 (h (µ−1) l )p+1 ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )ql(s)ds (27) A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 16 of 29 + Nµ−1∑ l=0 µ−1∑ v=0 h (v) l ∫ 1 0 K(t (µ) n,i , t (v) l + sh (v) l )dse(t (v) l ) + Nµ−1∑ l=0 µ−1∑ v=0 r−1∑ k=0 h (v) l ∫ 1 0 K(t (µ) n,i , t (v) l + sh (v) l )h (v) l αk(s)dse ′(v) l−k + Nµ−1∑ l=0 µ−1∑ v=0 h (v) l ∫ 1 0 K(t (µ) n,i , t (v) l + sh (v) l )h (v) l m∑ j=1 βj(s)dsE (v) l,j + n−1∑ l=r h (µ) l ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )dse(t (µ) l ) + n−1∑ l=r r−1∑ k=0 (h (µ) l )2 ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )αk(s)dse ′(µ) l−k + n−1∑ l=r m∑ j=1 (h (µ) l )2 ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )βj(s)dsE (µ) l,j +h (µ) n ∫ si 0 K(t (µ) n,i , t (µ) n + sh(µ)n )dse(t(µ)n ) +(h (µ) n )2 ∫ si 0 K(t (µ) n,i , t (µ) n + sh(µ)n ) r−1∑ k=0 αk(s)dse ′(µ) n−k +(h (µ) n )2 ∫ si 0 K(t (µ) n,i , t (µ) n + sh(µ)n ) m∑ j=1 βj(s)dsE (µ) n,j + Nµ−1∑ l=0 µ−2∑ v=0 h (v) l ∫ 1 0 K̂(t (µ) n,i , t (v) l + sh (v) l )dse(t (v) l ) + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )2 r−1∑ k=0 ∫ 1 0 K̂(t (µ) n,i , t (v) l + sh (v) l )αk(s)dse ′(v) l−k + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )2 m∑ j=1 ∫ 1 0 K̂(t (µ) n,i , t (v) l + sh (v) l )βj(s)dsE (v) l,j + n−1∑ l=r h (µ−1) l ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )dse(t (µ−1) l ) + n−1∑ l=r (h (µ−1) l )2 r−1∑ k=0 ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )αk(s)dse ′(µ−1) l−k + n−1∑ l=r (h (µ−1) l )2 m∑ j=1 ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )βj(s)dsE (µ−1) l,j +h (µ−1) n ∫ si 0 K̂(t (µ) n,i , t (µ−1) n + sh(µ−1) n )dse(t(µ−1) n ) +(h (µ−1) n )2 r−1∑ k=0 ∫ si 0 K̂(t (µ) n,i , t (µ−1) n + sh(µ−1) n )αk(s)dse ′(µ−1) n−k +(h (µ−1) n )2 m∑ j=1 ∫ si 0 K̂(t (µ) n,i , t (µ−1) n + sh(µ−1) n )βj(s)dsE (µ−1) n,j . A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 17 of 29 Now, we define the following matrices and vectors: εεε (µ) n = ( e ′(µ) n , ..., e ′(µ) n−r+1 )T , EEE (µ) n = ( E (µ) n,1 , ..., E (µ) n,m )T , DDD (.),l n =  (∫ 1 0 K(t (µ) n,i , t (.) l + sh (.) l )β1(s)ds, ..., ∫ 1 0 K(t (µ) n,i , t (.) l + sh (.) l )βm(s)ds ) , l = r − 1, . . . , n− 1, (∫ si 0 K(t (µ) n,i , t (.) l + sh (.) l )β1(s)ds, ..., ∫ si 0 K(t (µ) n,i , t (.) l + sh (.) l )βm(s)ds ) , l = n. D̃DD (.),l n =  (∫ 1 0 K̂(t (µ) n,i , t (.) l + sh (.) l )β1(s)ds, ..., ∫ 1 0 K̂(t (µ) n,i , t (.) l + sh (.) l )βm(s)ds ) , l = r − 1, . . . , n− 1, (∫ si 0 K̂(t (µ) n,i , t (.) l + sh (.) l )β1(s)ds, ..., ∫ si 0 K̂(t (µ) n,i , t (.) l + sh (.) l )βm(s)ds ) , l = n. GGG (.),l n =  (∫ 1 0 K(t (µ) n,i , t (.) l + sh (.) l )α0(s)ds, ..., ∫ 1 0 K(t (µ) n,i , t (.) l + sh (.) l )αr−1(s)ds ) , l = r − 1, . . . , n− 1, (∫ si 0 K(t (µ) n,i , t (.) l + sh (µ) l )α0(s)ds, ..., ∫ si 0 K(t (µ) n,i , t (.) l + sh (.) l )αr−1(s)ds ) , l = n. G̃GG (.),l n =  (∫ 1 0 K̂(t (µ) n,i , t (.) l + sh (.) l )α0(s)ds, ..., ∫ 1 0 K̂(t (µ) n,i , t (.) l + sh (.) l )αr−1(s)ds ) , l = r − 1, . . . , n− 1, (∫ si 0 K̂(t (µ) n,i , t (.) l + sh (µ) l )α0(s)ds, ..., ∫ si 0 K̂(t (µ) n,i , t (.) l + sh (.) l )αr−1(s)ds ) , l = n. HHH (.) n,κ,ϖ = ( h(.)κ (∫ 1 0 K(t (µ) n,i , t (.) κ + sh(.)κ )ds, ..., h(.)ϖ ∫ 1 0 K(t (µ) n,i , t (.) ϖ + sh(.)ϖ )ds ) , H̃HH (.) n,κ,ϖ = ( h(.)κ (∫ 1 0 K̂(t (µ) n,i , t (.) κ + sh(.)κ )ds, ..., h(.)ϖ ∫ 1 0 K̂(t (µ) n,i , t (.) ϖ + sh(.)ϖ )ds ) , LLL (.) 1,n = diag (∫ s1 0 K(t (µ) n,1, t (.) n + sh(.)n )ds, ..., ∫ sm 0 K(t(µ)n,m, t(.)n + sh(.)n )ds ) , L̃LL (.) 1,n = diag (∫ s1 0 K̂(t (µ) n,1, t (.) n + sh(.)n )ds, ..., ∫ sm 0 K̂(t(µ)n,m, t(.)n + sh(.)n )ds ) , ααα = ( α0(ci), ..., αr−1(ci) ) , βββ = ( β1(ci), ..., βm(ci) ) , ΨΨΨ (µ,ν),(p+1,p+2) m,n,l = (Ψ (µ,ν),(p+1,p+2) m,n,l )i. (28) A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 18 of 29 where for i = 1, · · · ,m, we have (Ψ (µ,ν),(p+1,p+2) m,n,l )i = (h (µ) n )p+1c1(t (µ) n,i )R (µ) m+1,n(si) +(h (µ−1) n )p+1c2(t (µ) n,i )R (µ−1) m+1,n(si) + Nµ−1∑ l=0 µ−1∑ v=0 (h (v) l )p+2 ∫ 1 0 K(t (µ) n,i , t (v) l + sh (v) l )R (v) m+1,l(s)ds + n−1∑ l=r (h (µ) l )p+2 ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l )R (µ) m+1,l(s)ds +(h (µ) n )p+2 ∫ si 0 K(t (µ) n,i , t (µ) n + sh(µ)n )R (µ) m+1,n(s)ds + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )p+2 ∫ 1 0 K̂(t (µ) n,i , t (v) l + sh (v) l )R (v) m+1,l(s)ds + n−1∑ l=r (h (µ−1) l )p+2 ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l ) +R (µ−1) m+1,l(s)ds +(h (µ−1) n )p+2 ∫ si 0 K̂(t (µ) n,i , t (µ−1) n + sh(µ−1) n )R (µ−1) m+1,n(s)ds + r−1∑ l=1 (h (µ) l )p+1 ∫ 1 0 K(t (µ) n,i , t (µ) l + sh (µ) l ql(s)ds + r−1∑ l=1 (h (µ−1) l )p+1 ∫ 1 0 K̂(t (µ) n,i , t (µ−1) l + sh (µ−1) l )ql(s)ds. (29) By substituting these matrices and vectors into the equation (28), the matrix form of A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 19 of 29 the equation (27), can be written as follows: [ III − h (µ) n ( CCC (µ) 1,nβββ + h (µ) n DDD (µ),n n ) −h (µ) n ( CCC (µ) 1,nααα+ h (µ) n GGG (µ),n n )]EEE(µ) n εεε (µ) n  = h (µ) n [ CCC (µ) 2,nβββ 000 ] [EEE(µ−1) n εεε (µ−1) n ] + Nµ−1∑ l=0 µ−1∑ v=0 (h (v) l )2 [ DDD (ν),l n 000 ] [EEE(v) l εεε (v) l ] + n−1∑ l=r (h (µ) l )2 [ DDD (µ),l n 000 ] [EEE(µ) l εεε (µ) l ] + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )2 [ D̃DD (ν),l n 000 ] [EEE(v) l εεε (v) l ] + n−1∑ l=r (h (µ−1) l )2 [ D̃DD (µ−1),l n 000 ] [EEE(µ−1) l εεε (µ−1) l ] + (h(µ−1) n )2 [ D̃DD (µ−1),n n 000 ] [EEE(µ−1) n εεε (µ−1) n ] + Nµ−1∑ l=0 µ−1∑ v=0 (h (v) l )2 [ 000 GGG (v),l n ] [EEE(v) l ε (v) l ] + n−1∑ l=r (h (µ) l )2 [ 000 GGG (µ),l n ] [EEE(µ) l ε (µ) l ] + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )2 [ 000 G̃GG (v),l n ] [EEE(v) l ε (v) l ] + n−1∑ l=r (h (µ−1) l )2 [ 000 G̃GG (µ−1),l n ] [EEE(µ−1) l ε (µ−1) l ] +(h (µ−1) n )2 [ 000 G̃GG (µ−1),n n ] [EEE(µ−1) n ε (µ−1) n ] + h (µ) n [ 000 CCC (µ) 2,nα̃αα ] [EEE(µ−1) n ε (µ−1) n ] +FFF (v) +ΨΨΨ (µ,ν),(p+1,p+2) m,n,l , (30) where FFF (v) = µ−1∑ v=0 HHH (v) n,0,Nµ−1ϵϵϵ (v) 0,Nµ−1 +HHH (µ) n,r,n−1ϵϵϵ (µ) r,n−1 + µ−2∑ v=0 H̃HH (v) n,0,Nµ−1ϵϵϵ (v) 0,Nµ−1 + h(µ)n LLL (µ) 1,ne(t (µ) n ) + H̃HH (µ−1) n,r,n−1ϵϵϵ (µ−1) r,n−1 +h (µ−1) n L̃LL (µ−1) 1,n e(t (µ−1) n ) +CCC (µ) 1,ne(t (µ) n ) +CCC (µ) 2,ne(t (µ−1) n ). (31) Also, equation (26) with z = 1, leads to εεε (.) l = PPPεεε (.) l−1 +QQQEEE (.) l−1 +O(h(.))p, (32) where PPP = [ 000r−1,1 Ir−1 Pr−1(1) Pr−2(1), · · · , P0(1), ] , QQQ =  000r−1,m Q1(1) · · · Qm(1)  . A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 20 of 29 The solution of the difference equation (32) is εεε (.) l = PPP l−r+1εεε (.) r−1 + l−1∑ j=r−1 PPP l−j−1QQQE (.) j +O(h(.))p. (33) Now, from equations (30) and (32), we get[ III − h (µ) n ( CCC (µ) 1,nβββ + h (µ) n DDD (µ),n n ) −h (µ) n ( CCC (µ) 1,nααα+ h (µ) n GGG (µ),n n ) 000 III ][ EEE (µ) n εεε (µ) n ] = [ (h (µ) n−1) 2DDD (µ),n−1 n 000 QQQ PPP ][ EEE (µ) n−1 εεε (µ) n−1 ] +h (µ) n [ CCC (µ) 2,nβββ 000 000 000 ][ EEE (µ−1) n εεε (µ−1) n ] + Nµ−1∑ l=0 µ−1∑ v=0 (h (v) l )2 [ DDD (ν),l n 000 000 000 ][ EEE (v) l εεε (v) l ] + n−2∑ l=r (h (µ) l )2 [ DDD (µ),l n 000 000 000 ][ EEE (µ) l εεε (µ) l ] + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )2 [ D̃DD (ν),l n 000 000 000 ][ EEE (v) l εεε (v) l ] + n−1∑ l=r (h (µ−1) l )2 [ D̃DD (µ−1),l n 000 000 000 ][ EEE (µ−1) l εεε (µ−1) l ] + (h(µ−1) n )2 [ D̃DD (µ−1),n n 000 000 000 ][ EEE (µ−1) n εεε (µ−1) n ] + Nµ−1∑ l=0 µ−1∑ v=0 (h (v) l )2 [ 000 GGG (v),l n 000 000 ][ EEE (v) l ε (v) l ] + n−1∑ l=r (h (µ) l )2 [ 000 GGG (µ),l n 000 000 ][ EEE (µ) l ε (µ) l ] + Nµ−1∑ l=0 µ−2∑ v=0 (h (v) l )2 [ 000 G̃GG (v),l n 000 000 ][ EEE (v) l ε (v) l ] + n−1∑ l=r (h (µ−1) l )2 [ 000 G̃GG (µ−1),l n 000 000 ][ EEE (µ−1) l ε (µ−1) l ] +(h (µ−1) n )2 [ 000 G̃GG (µ−1),n n 000 000 ][ EEE (µ−1) n ε (µ−1) n ] + h (µ) n [ 000 CCC (µ) 2,nα̃αα 000 000 ][ EEE (µ−1) n ε (µ−1) n ] + FFF (v) 000 + [ ΨΨΨ (µ,ν),(p+1,p+2) m,n,l O(h (µ) n )p ] , (34) where expression (25), with z = 1, leads to e(t (.) n ) = e(t (.) n−1 + h (.) n ) = e(t (.) n−1) + h (.) n−1 r−1∑ k=0 αk(1)e ′(.) n−1−k + h (.) n−1 m∑ j=1 βj(1)E (.) n−1,j + (h (.) n−1) p+1R (.) m+r,n−1(1). A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 21 of 29 Note that, in the equation (34), the matrix[ III − h (µ) n ( CCC (µ) 1,nβββ + h (µ) n DDD (µ),n n ) −h (µ) n ( CCC (µ) 1,nααα+ h (µ) n GGG (µ),n n ) 000 III ] , µ = 0, 1, ..., coincides with Theorem 4.5.3 in [1] and Theorem 5.1 in [36] and then this matrix is invertible, so it’s inverse is uniformly bounded. Now, with assumption h = max l,ν h (ν) l , we have ∥ E(1) ∥∞≤ C1h p, where C1 is a constant and E(1) = [ EEE (µ) n εεε (µ) n ] . Hence, this proves Theorem 3, in same manner in Theorems 3.2.3 and 4.5.2, in [1]. Remark 4. We conclude this section with a comment regarding the extension of the results of Theorem 3 to the non-linear equation (1). Under the assumption of the existence of a (unique) solution x(t) on J , the non-linear analogue of the error equation (1) is e′(t) = c1(t)e(t) + c2(t)e(τ(t)) + δ(t) + ∫ t t0 ( K(t, s, x(s))−K(t, s, w(s)) ) ds+ ∫ τ(t) t0 ( K̂(t, s, x(s))− K̂(t, s, w(s)) ) ds. (35) If the partial derivatives ∂K ∂x and ∂K̂ ∂x are continuous and bounded on the domain of their own definition. Assuring the existence of a unique collocation solution w, then (35) may again be written in the form (17). The roles of K and K̂ are now assumed by k(t, s) = ∂K(t, s, Z(s)) ∂x , and k̂(t, s) = ∂K̂(t, s, Z(s)) ∂x , where Z(s) = θx(s) + (1 − θ)w. Hence, the above proof is easily adapted to deal with the non-linear case (1), and so the convergence results of Theorem 3 remain valid for non-linear equations. 4. Numerical examples Here, we examine two numerical instances to demonstrate the effectiveness of the suggested method. All calculations were executed using Mathematica® software, Version 11.1. We choose s1 = 0.8, s2 = 1 and r = 2. Since s2 = 1, then we have ρ(PPP ) = 0 < 1. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 22 of 29 Throughout the subsequent part of this section, all numerical experiments employ T = 1, and the initial values are derived from well-established exact solutions. The analysis of the numerical results presented in the Tables reveals that the multi-step method demonstrates greater accuracy compared to the one-step method employed in [1]. In Tables 1,2, we report the maximum of the absolute errors at the grid points for m = 2 and r = 2. Also, we calculate the order of convergence by p = log2 ( ∥eN∥∞ ∥e2N∥∞ ) , and report it in Table 3. From Theorem 3, we know that for m, r = 2, the order of convergence is equal to p = m+ r = 4 while for one-step collocation methods, the order of convergence is equal to p = 2. Additionally, in Figures 1 and 2, we graph the convergence order for both the one-step and multi-step schemes across various values of N . For the one-step method, order of convergence tends to p = 2 and for multi-step schemes, that tends to p = 4. The observed order of convergence aligns with the theoretical findings outlined in Theorem 3. Obviously, noting Tables 1, 2 and Figures 1, 2, we see that using the multi-step collocation methods can get higher convergence orders than the classical one-step collocation methods when the same number of collocation parameters is used. Example 1. Examine VIDEs with non-vanishing delay x′(t) = −tx(t)− (t+ 1)x(12 t) + f(t) + ∫ t 1 4 sin(s− t)x2(s)ds+ ∫ 1 2 t 1 4 t cos(s)x2(s)ds, t ∈ [14 , 1], x(t) = e2t, t ∈ [18 , 1 4 ], and f so that the exact solution is x(t) = e2t. We use FindRoot command in Mathematica software to solve non-linear algebraic equations associated with the nonlinear test problem 1. If you specify only one starting value, FindRoot searches for a solution using Newton methods. If FindRoot does not succeed in finding a solution to the accuracy you specify within MaxIterations steps, it returns the most recent approximation to a solution that it found. You can then apply FindRoot again, with this approximation as a starting point. If we do not choose the starting point correctly, the software warns when running and we can change the starting point. Example 2. Examine VIDEs with non-vanishing delay x′(t) = (t2 + 2)x(12 − t) + t+ ∫ t− 1 2 0 (2s+ 3t+ 1)x(s)ds, t ∈ [0, 1], x(t) = 1, t ∈ [−1 2 , 0], A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 23 of 29 and the exact solution is x(t) =  1 + 7t 4 − t2 4 + 5t3 3 , 0 < t ≤ 1 2 , 6847 4608 − 59t 128 + 433t2 128 − 299t3 144 + 109t4 32 − 11t5 8 + 43t6 72 , 1 2 < t ≤ 1. Resulting from the delay functionτ(t) = t− 1 2 , we have ςµ = µ 1 2 , µ = 0, 1. Also lim t→0− x′(t) ̸= lim t→0+ x′(t), lim t→ 1 2 − x′(t) = lim t→ 1 2 + x′(t). Table 1: L∞ errors and CPU time based on seconds for m, r = 2 in Example 1. (s1, s2) = (0.8, 1) (s1, s2) = (0.8, 1) N One-step method [1] CPU time(sec) Multi-step method CPU time(sec) 4 2.74× 10−2 0.391 2.56× 10−5 1.64 8 6.94× 10−3 0.563 2.14× 10−6 2.14 16 1.74× 10−3 1.29 1.47× 10−7 9.46 32 4.37× 10−4 4.12 9.64× 10−9 43.6 Table 2: L∞ errors and CPU time based on seconds for m, r = 2 in Example 2. (s1, s2) = (0.8, 1) (s1, s2) = (0.8, 1) N One-step method [1] CPU time(sec) Multi-step method CPU time(sec) 4 5.14× 10−2 0.469 2.06× 10−5 0.672 8 1.24× 10−3 0.562 1.74× 10−6 2.70 16 3.06× 10−3 1.23 1.20× 10−7 7.79 32 7.58× 10−4 4.36 7.86× 10−9 29.9 Table 3: Order of convergence for m = r = 2 in Examples 1 and 2. Order of convergence Order of convergence N One-step for Ex. 1 Multi-step for Ex. 1 One-step for 2 Multi-step for Ex. 2 8 1.884 3.545 2.146 3.559 16 1.993 3.838 2.024 3.855 32 1.996 3.929 2.012 3.941 64 1.999 3.989 2.000 3.991 A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 24 of 29 10 20 30 40 50 60 1 2 3 4 5 6 N O rd e r o f c o n v e rg e n c e one-step Multi-step Figure 1: Plot of order of convergence for the one-step and multi-step scheme in Example 1. 10 20 30 40 50 60 1 2 3 4 5 6 N O rd e r o f c o n v e rg e n c e one-step Multi-step Figure 2: Plot of order of convergence for the one-step and multi-step scheme in Example 2. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 25 of 29 In this section, by considering m = r = 2, the value of n begins at r − 1 = 1 for the approximate solutions (12). According to the proposed numerical methods, we need starting values x ′(µ) 0 = x′(t (µ) 0 ), x ′(µ) 1 = x′(t (µ) 1 ) = x′(t (µ) 0 + h (µ) 0 ), x (µ) 1 = x(t (µ) 1 ) = x(t (µ) 0 + h (µ) 0 ) and w(t (µ) 0 + sh (µ) 0 ) for µ = 0, 1. Every value in the above list was derived from the known exact solutions. We can utilize approximate solutions for the starting values because, as you are aware, we do not have the exact solutions for real problems. The following Remark is taken into consideration for this purpose: Remark 5. We can examine the impact of numerical approximations of the initial values using a traditional one-step approach. We consider the polynomial approximation for the equation (1) as follows, based on [1]: w′(t(µ)n + zh(µ)n ) = h(µ)n m∑ j=1 Lj(z)W (µ) n,j , n = 0, . . . , N − 1, (36) and w(t(µ)n + zh(µ)n ) = x(µ)n + h(µ)n m∑ j=1 βj(z)W (µ) n,j , n = 0, . . . , N − 1, (37) where βj(z) = ∫ z 0 Lj(s)ds and Lj(j = 1, . . . ,m) represent the Lagrange canonical polyno- mials. Taking into account Theorem 4.5.2 [1], we have the collocation error related to the collocation solution ∥x(ν) − w(ν)∥∞ ≤ Cνh m, ν = 0, 1, (38) where the constants Cν are independent of h. For µ = 0, 1, the approximation of the starting values x ′(µ) 1 = x′(t (µ) 0 + h (µ) 0 ) and x (µ) 1 = x(t (µ) 0 + h (µ) 0 ) can be obtained from (36) and (37) by setting z = 1. It is also possible to obtain w(t (µ) 0 + sh (µ) 0 ) from (37). We used the multi-step collocation method with m = 2 and r = 2 and looked at two numerical cases. Next, using the approximation technique with the order m+r = 4, we must determine the initial values based on Theorem 3. The starting values are obtained using a numerical technique based on the suggested classical one step methods (37) with the order 4 (see (38)) and collocation parameters si = 1 5−i , (i = 1, . . . , 4). Table 4 reports the maximum errors for various values of N . Table 4: L∞ errors with the approximate starting values and m = r = 2. N Multi-step method for Example 1 Multi-step method for Example 2 4 9.12× 10−5 9.01× 10−5 8 8.17× 10−6 8.56× 10−6 16 7.66× 10−7 8.18× 10−7 32 1.30× 10−8 1.81× 10−8 We also solve Example 1 and 2 with r = 3, m = 2 and report the results in the Table 5. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 26 of 29 Table 5: L∞ errors and CPU time based on seconds for m = 2 and r = 3 in Examples 1 and 2. (s1, s2) = (0.8, 1) (s1, s2) = (0.8, 1) N Multi-step for Example 1 CPU time(sec) Multi-step for Example 2 CPU time(sec) 4 1.24× 10−6 1.90 1.27× 10−6 1.26 8 1.09× 10−7 2.71 9.96× 10−8 2.93 16 4.11× 10−9 10.2 4.04× 10−9 9.96 32 1.38× 10−10 77.9 1.41× 10−10 40.6 5. Conclusion We have shown that the multi-step collocation method offers a reliable and precise numerical technique for estimating solutions to VIDEs with non-zero delay. Finally, to substantiate the theoretical predictions in a practical context, we examined some test problems. These tests served as evidence of the consistency and agreement between the numerical and theoretical analysis. We have considered our proposed methods based on the smoothness of the given function in the Theorems 1 and 2. As future work, for the case the solution derivatives are unbounded especially when either the data is non smooth may be able to use adaptive generated meshes in [37–39]. Also, we will study the proposed multi-step method to solve delay integro-differential-algebraic equations in the following form: A(t)X ′(t) = F (t) +B(t)X(t) + ∫ t 0 K(t, s,X(s))ds+ ∫ τ(t) 0 K̂(t, s,X(s))ds, t ∈ J, subject to detA(t) = 0, ∀t ∈ J. Unlike integro-differential algebraic equations (IDAEs), the solutions of delay integro-differential algebraic equations (DIDAEs) can be included primary discontinuity points. We will try to investigate the structure of the solution of DIDAEs from numerical and theoretical point of view. References [1] H. Brunner. Collocation methods for Volterra integral and related functional differ- ential equations. Cambridge University Press, Cambridge, 2004. [2] F. Brauer. On a nonlinear integral equation for population growth problems. SIAM Journal on Mathematical Analysis, 6(2):312–317, 1975. [3] F. Brauer and C. Castillo-Chavez. Mathematical models in population biology and epidemiology. Springer, New York, 2001. [4] B. Chen, J. Hu, and B. K. Ghosh. Finite-time tracking control of heterogeneous multi- AUV systems with partial measurements and intermittent communication. Science China Information Sciences, 67(5):152202, 2024. [5] Y. Kai, J. Ji, and Z. Yin. Study of the generalization of regularized long-wave equa- tion. Nonlinear Dynamics, 107(3):2745–2752, 2022. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 27 of 29 [6] L. Liu, S. Zhang, L. Zhang, G. Pan, and J. Yu. Multi-UUV maneuvering counter- game for dynamic target scenario based on fractional-order recurrent neural network. IEEE Transactions on Cybernetics, 53(6):4015–4028, 2023. [7] G. A. Bocharov and F. A. Rihan. Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics, 125(1- 2):183–199, 2000. [8] K. Maleknejad, F. Mirzaee, and S. Abbasbandy. Solving linear integro-differential equations system by using rationalized Haar functions method. Applied Mathematics and Computation, 155(2):317–328, 2004. [9] H. Brunner. Volterra integral equations: an introduction to theory and applications. Cambridge University Press, Cambridge, 2017. [10] H. Brunner, A. Makroglou, and R. K. Miller. Mixed interpolation collocation methods for first and second order Volterra integro-differential equations with periodic solution. Applied Numerical Mathematics, 23(4):381–402, 1997. [11] M. R. Crisci, E. Russo, and A. Vecchio. Stability of collocation methods for Volterra integro-differential equations. Journal of Integral Equations and Applica- tions, 4(4):491–507, 1992. [12] T. Lin, Y. Lin, M. Rao, and S. Zhang. Petrov–Galerkin methods for linear Volterra integro-differential equations. SIAM Journal on Numerical Analysis, 38(3):937–963, 2000. [13] Y. Jafarzadeh and B. Keramati. Numerical method for a system of integro-differential equations by Lagrange interpolation. Asian-European Journal of Mathematics, 9(3):1650058, 2016. [14] A. Cardone and D. Conte. Multistep collocation methods for Volterra integro- differential equations. Applied Mathematics and Computation, 221:770–785, 2013. [15] A. Akyüz and M. Sezer. Chebyshev polynomial solutions of systems of higher-order linear Fredholm-Volterra integro-differential equations. Journal of the Franklin Insti- tute, 342(6):688–701, 2005. [16] A. Karamete and M. Sezer. A Taylor collocation method for the solution of lin- ear integro-differential equations. International Journal of Computer Mathematics, 79(9):987–1000, 2002. [17] J. T. Katsikadelis. Numerical solution of integrodifferential equations with convolu- tion integrals. Archive of Applied Mechanics, 89(10):2019–2036, 2019. [18] K. Kumar, P. C. Podila, P. Das, and H. Ramos. A graded mesh refinement approach for boundary layer originated singularly perturbed time-delayed parabolic convection diffusion problems. Mathematical Methods in the Applied Sciences, 44(16):12332– 12350, 2021. [19] S. Saini, P. Das, and S. Kumar. Parameter uniform higher order numerical treat- ment for singularly perturbed Robin type parabolic reaction diffusion multiple scale problems with large delay in time. Applied Numerical Mathematics, 196:1–21, 2024. [20] P. Das. An a posteriori based convergence analysis for a nonlinear singularly perturbed system of delay differential equations on an adaptive mesh. Numerical Algorithms, 81(2):465–487, 2019. A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 28 of 29 [21] I. Ali, H. Brunner, and T. Tang. Spectral methods for pantograph-type differential and integral equations with multiple delays. Frontiers of Mathematics in China, 4(1):49–61, 2009. [22] H. Brunner. Collocation and continuous implicit Runge-Kutta methods for a class of delay Volterra integral equations. Journal of Computational and Applied Mathemat- ics, 53(1):61–72, 1994. [23] H. Brunner and Y. Yatsenko. Spline collocation methods for nonlinear Volterra in- tegral equations with unknown delay. Journal of Computational and Applied Mathe- matics, 71(1):67–81, 1996. [24] M. V. Bulatov, M. N. Machkhina, and V. N. Phat. Existence and uniqueness of solutions to nonlinear integral-algebraic equations with variable limits of integrations. Communications on Applied Nonlinear Analysis, 21(1):65–76, 2014. [25] F. Caliò, E. Marchetti, and R. Pavani. About the deficient spline collocation method for particular differential and integral equations with delay. Rendiconti del Seminario Matematico, 61(3):287–300, 2003. [26] F. Caliò, E. Marchetti, R. Pavani, and G. Micula. About some Volterra problems solved by a particular spline collocation. Studia Universitatis Babeş-Bolyai Mathe- matica, 48(2):45–52, 2003. [27] M. Khasi, F. Ghoreishi, and M. Hadizadeh. Numerical analysis of a high order method for state-dependent delay integral equations. Numerical Algorithms, 66(1):177–201, 2014. [28] S. Kumar, S. Kumar, and Sumita. A priori and a posteriori error estimation for singularly perturbed delay integro-differential equations. Numerical Algorithms, 95(4):1561–1582, 2024. [29] S. Wu and S. Gan. Errors of linear multistep methods for singularly perturbed Volterra delay integro-differential equations. Mathematics and Computers in Simu- lation, 79(10):3148–3159, 2009. [30] N. A. Baharum, Z. A. Majid, N. Senu, and H. Rosali. Numerical approach for delay Volterra integro-differential equation. Sains Malaysiana, 51(12):4125–4144, 2022. [31] B. Gürbüz. A numerical scheme for the solution of neutral integro-differential equa- tions including variable delay. Mathematical Sciences, 16:13–21, 2022. [32] Q. Guan, R. Zhang, and Y. Zou. Analysis of collocation methods for nonstandard Volterra integral equations. IMA Journal of Numerical Analysis, 32(4):1755–1785, 2012. [33] M. Chandru, P. Das, and H. Ramos. Numerical treatment of two-parameter singularly perturbed parabolic convection diffusion problems with non-smooth data. Mathemat- ical Methods in the Applied Sciences, 41(14):5359–5387, 2018. [34] P. Das and S. Natesan. Optimal error estimate using mesh equidistribution tech- nique for singularly perturbed system of reaction–diffusion boundary-value problems. Applied Mathematics and Computation, 249:265–277, 2014. [35] S. Saini, P. Das, and S. Kumar. Computational cost reduction for coupled system of multiple scale reaction diffusion problems with mixed type boundary conditions having boundary layers. Revista de la Real Academia de Ciencias Exactas, F́ısicas y A. Ali Eashel, S. Pishbin, P. Darania / Eur. J. Pure Appl. Math, 18 (2) (2025), 5766 29 of 29 Naturales - Serie A: Matemáticas, 117(2):66, 2023. [36] P. Darania and S. Pishbin. Multistep collocation methods for integral-algebraic equa- tions with non-vanishing delays. Mathematics and Computers in Simulation, 205:33– 61, 2023. [37] P. Das, S. Rana, and H. Ramos. On the approximate solutions of a class of fractional order nonlinear Volterra integro-differential initial value problems and boundary value problems of first kind and their convergence analysis. Journal of Computational and Applied Mathematics, 404:113116, 2022. [38] S. Kumar, Sumita, and J. Vigo-Aguiar. A high order convergent numerical method for singularly perturbed time dependent problems using mesh equidistribution. Math- ematics and Computers in Simulation, 199:287–306, 2022. [39] S. Santra, J. Mohapatra, P. Das, and D. Choudhuri. Higher order approximations for fractional order integro-parabolic partial differential equations on an adaptive mesh with error analysis. Computers & Mathematics with Applications, 150:87–101, 2023.