IHJPAS. 36 (3) 2023 382 This work is licensed under a Creative Commons Attribution 4.0 International License *Corresponding Author: ghada.h.i@ihcoedu.uobaghdad.edu.iq Abstract This work describes two efficient and useful methods for solving fractional pantograph delay equations (FPDEs) with initial and boundary conditions. These two methods depend mainly on orthogonal polynomials, which are the method of the operational matrix of the fractional derivative that depends on Bernstein polynomials and the operational matrix of the fractional derivative with Shifted Legendre polynomials. The basic procedure of this method is to convert the pantograph delay equation to a system of linear equations, and by using, the operational matrices, we get rid of the integration and differentiation operations, which makes solving the problem easier. The concept of Caputo has been used to describe fractional derivatives. Finally, some numerical examples are identified to show the utility and capability of the two proposed approaches. The Mathematica® 12 program has been relied upon in the calculations. Keywords: Bernstein polynomials, Legendre polynomials, Pantograph Delay Equations. 1. Introduction Fractional calculus has recently been increasingly used due to its wide applications in real-world problems. It has been used as a powerful tool in various fields such as chemistry, physics, engineering, and applied sciences, where it can be seen in thermal modeling [1], networks [2], optics [3], optimal control [4], elasticity [5], fluid mechanics [6], and many other applications. The fractional nature of this class of problems makes the solution very difficult. As a result, many researchers and authors try to generalize and develop the current methods to simply apply them numerically or analytically and find approximate solutions to them. These methods include the Collocation method [7, 8], Adomian’s decomposition method [9–11], homotopy perturbation method [12, 13], and other methods [14–16]. Fractional delay differential equations (FDDEs) are used in a variety of technical systems, including characterizing propagation, transport development, or population gestures [17, 18]. The pantograph equation is one of the most important types of delay differential equations, and it is used to describe a wide range of phenomena. Various applications of these equations in practical doi.org/10.30526/36.3.3076 Article history: Received 16 October 2022, Accepted 9 January 2023, Published in July 2023. Ibn Al-Haitham Journal for Pure and Applied Sciences Journal homepage: jih.uobaghdad.edu.iq Fractional Pantograph Delay Equations Solving by the Meshless Methods Shefaa M. N. Jasim Department of Mathematics, College of Education for Pure Science Ibn Al-Haytham, University of Baghdad, Baghdad, Iraq. Shefaa.Mohammed1203a@ihcoedu.uobaghdad.edu.i q Ghada H. Ibraheem* Department of Mathematics, College of Education for Pure Science Ibn Al-Haytham, University of Baghdad, Baghdad, Iraq. ghada.h.i@ihcoedu.uobaghdad.edu.iq https://creativecommons.org/licenses/by/4.0/ mailto:ghada.h.i@ihcoedu.uobaghdad.edu.iq http://en.uobaghdad.edu.iq/?page_id=15060 http://en.uobaghdad.edu.iq/?page_id=15060 mailto:Shefaa.Mohammed1203a@ihcoedu.uobaghdad.edu.iq mailto:Shefaa.Mohammed1203a@ihcoedu.uobaghdad.edu.iq mailto:Shefaa.Mohammed1203a@ihcoedu.uobaghdad.edu.iq http://en.uobaghdad.edu.iq/?page_id=15060 http://en.uobaghdad.edu.iq/?page_id=15060 mailto:ghada.h.i@ihcoedu.uobaghdad.edu.iq mailto:ghada.h.i@ihcoedu.uobaghdad.edu.iq IHJPAS. 36 (3) 2023 383 disciplines such as biology, physics, economics, and electrodynamics have been researched by many scholars [19–21]. Researchers have been interested in finding strategies to solve these kinds of problems and have presented a number of papers on the subject, where Rabiei and Ordok [22] introduce fractional-order Boubaker polynomials related to the Boubaker polynomials to achieve the numerical result for pantograph differential equations of fractional order in any arbitrary interval. Yuttanan et al. [23] presented a new class of functions called fractional-order generalized Taylor wavelets (FOGTW) for the nonlinear fractional delay and nonlinear fractional pantograph differential equations. In [24], a method utilizing the Chebyshev cardinal functions (CCFs) is formulated to find an accurate result for fractional pantograph delay differential equations; see more works [25–27]. Operational matrices are one of the most widely used numerical methods for solving fractional differential equations (FDEs) because they are based on polynomials, which are one of the simplest functions in terms of structure. As a result, the benefit of this method is that it transforms the problem into a set of algebraic equations. Using the operational matrix, we get rid of the derivations and integrals, making the solution easier and simpler. Studies that relied on operational matrices methods can be seen in [28–32]. The main objective of this work is to apply the operational matrix method for derivatives that depend on Bernstein polynomials, as well as the operational matrix method, which depends on shifted Legendre polynomials, for solving fractional pantograph delay equations (FPDEs) with initial and boundary conditions. The rest of this paper is organized as follows: In Section 2, we provide some basic definitions of fractional calculus. In Section 3, Bernstein polynomials and their operational matrices (BOM) are introduced. Also, in Section 4, Shifted Legendre polynomials and their operational matrices (LOM) are considered. In Section 5, some numerical examples of pantograph equations with prime and boundary conditions will be solved by the proposed methods. Finally, we discuss our findings and conclusions. 2. Basic Definitions The Riemann-Liouville and Caputo operators [33] are the two most commonly used definitions in fractional calculus. This section will give the definitions and some of their properties. Definition 1: The Riemann–Liouville fractional integration of order α is defined as [33] 𝐼𝛼 𝑦(𝑥) = { 1 𝛤(𝛼) ∫ 𝑥 0   𝑦(𝑠) (𝑥−𝑠)1−𝛼 𝑑𝑠, 𝛼 > 0, 𝑥 > 0, 𝑦(𝑥), 𝛼 = 0 . (1) For the Riemann–Liouville integral 1. 𝐼𝛼1 𝐼𝛼2 𝑦(𝑥) = 𝐼𝛼1+𝛼2 𝑦(𝑥), 2. 𝐼𝛼 (𝜆1𝑦(𝑥) + 𝜆2𝑔(𝑥)) = 𝜆1𝐼 𝛼 𝑦(𝑥) + 𝜆2𝐼 𝛼 𝑔(𝑥), 3. 𝐼𝛼 𝑥𝛽 = 𝛤(𝛽+1) 𝛤(𝛽+𝛼+1) 𝑥𝛼+𝛽 , 𝛽 > −1, where 𝛼1 , 𝛼2, 𝜆1 𝑎𝑛𝑑 𝜆2 are constants. Definition 2: Caputo’s fractional derivative operator of order 𝛼 is defined as [33] 𝐷𝛼 𝑦(𝑥) = 1 𝛤(𝑛−𝛼) ∫ 𝑥 0   𝑦(𝑛)(𝑠) (𝑥−𝑠)𝛼+1−𝑛 𝑑𝑠, 𝑛 − 1 < 𝛼 ≤ 𝑛, 𝑛 ∈ 𝑁. (2) IHJPAS. 36 (3) 2023 384 The characteristics of Caputo derivative are 1. 𝐷𝛼 𝐼𝛼 𝑦(𝑥) = 𝑦(𝑥), 2. 𝐼𝛼 𝐷𝛼 𝑦(𝑥) = 𝑦(𝑥) − ∑𝑛−1𝑖=0   𝑦 (𝑖)(0) 𝑥𝑖 𝑖! , 3. 𝐷𝛼 𝜆 = 0, Where 𝜆 is constant. In this study, the concept of Caputo was relied upon to define the fractional derivative. 3. Bernstein polynomials and their operational matrix [34,35] In this section, we will go over Bernstein's polynomial and some of its fundamental properties, as well as a description of an operational matrix with integer and fractional derivatives, and explain how to apply the approach (BOM) for solving the pantograph equation. 3. 1. Bernstein polynomials The Bernstein polynomials of degree 𝑛 on the interval [0,1] are defined by 𝐵𝑖,𝑛(𝑥) = (𝑛 𝑖 )𝑥 𝑖 (1 − 𝑥)𝑛−𝑖 , for 𝑖 = 0, 1, … , 𝑛, where (𝑛 𝑖 ) = 𝑛! 𝑖! (𝑛 − 𝑖)! . Bernstein basis polynomials in a linear combination 𝐵𝑛 (𝑥) = ∑ 𝑛 𝑖=0   𝑐𝑖 𝐵𝑖,𝑛(𝑥) , are called a Bernstein polynomial or polynomial in Bernstein form of degree , and 𝑐𝑖 are called Bernstein coefficients. Now we can approximate any polynomial of degree 𝑛 to the form of linear combination, as given below [34] 𝑦(𝑥) = ∑𝑛𝑖=1   𝑐𝑖 𝐵𝑖,𝑛 = 𝐶 𝑇 𝛷(𝑥), (3) where 𝐶𝑇 = [𝑐0, 𝑐1, … , 𝑐𝑛], 𝛷(𝑥) = [𝐵0,𝑛, 𝐵1,𝑛, … , 𝐵𝑛,𝑛], 𝑇 or in the form of a matrix resulting from multiplying a square matrix (𝑛 + 1) × (𝑛 + 1) and vector (𝑛 + 1) × 1 : 𝛷(𝑥) =𝐴𝑋, 𝐴 = [(−1)0(𝑛 0 ) (−1)1(𝑛 0 )(𝑛 − 0 1 ) … (−1)𝑛−0(𝑛 0 )(𝑛 − 0 𝑛 − 0 ) 0 (−1)0(𝑛 𝑖 ) … (−1)𝑛−𝑖 (𝑛 𝑖 )(𝑛 − 𝑖 𝑛 − 𝑖 ) ⋮ ⋮ ⋱ ⋮ 0 0 … (−1)0(𝑛 𝑛 ) ](𝑛+1)×(𝑛+1), 𝑋 = [1 𝑥 𝑥 2 ⋮ 𝑥𝑛 ] (𝑛+1)×1 . We note that the ‖𝐴‖ ≠ 0, which indicates that the matrix 𝐴 is invertible. On the other hand, when using the Bernstein polynomials in the least-squares approximation, the fact that they are not orthogonal becomes a disadvantage. According to [35], one method for direct least-squares approximation by polynomials in Bernstein form relies on the construction of the IHJPAS. 36 (3) 2023 385 basis {𝑑0 𝑛(𝑥), 𝑑1 𝑛(𝑥), … , 𝑑𝑛 𝑛(𝑥)} that is "dual" to the Bernstein basis of degree 𝑛 on to [0,1]. The attribute that distinguishes this dual basis is ∫ 1 0   𝑏𝑖 𝑛(𝑥)𝑑𝑗 𝑛(𝑥)𝑑𝑥 = {1 𝑓𝑜𝑟 𝑖 = 𝑗, 0 𝑓𝑜𝑟 𝑖 ≠ 𝑗, for 𝑖, 𝑗 = 0,1, . . . . . . . , 𝑛, where 𝑑𝑗 𝑛(𝑥) = ∑𝑛𝑘=0   𝜆𝑗𝑘 𝐵𝑘 𝑛(𝑥), 𝑗 = 0,1, … , 𝑛, 𝜆𝑗𝑘 = (−1)𝑗+𝑘 (𝑛 𝑗 )(𝑛 𝑘 ) ∑ 𝑚𝑖𝑛(𝑗,𝑘) 𝑖=0   (2𝑖 + 1)(𝑛 + 𝑖 + 1 𝑛 − 𝑗 )(𝑛 − 𝑖 𝑛 − 𝑗 )(𝑛 + 𝑖 + 1 𝑛 − 𝑘 )(𝑛 − 𝑖 𝑛 − 𝑘 ), (4) for 𝑗, 𝑘 = 0,1, . . . . , 𝑛. 3. 2 Operational Matrix of Fractional Derivative for Bernstein Polynomial (BOM) In the beginning, we will present the operational matrix 𝐷 of derivatives for Bernstein polynomials, which is a square matrix of degree (𝑛 + 1) × (𝑛 + 1), then the derivative of the vector 𝛷(𝑥) is ⅆ𝛷(𝑥) ⅆ𝑥 = 𝐷𝛷(𝑥), (5) where 𝐷 is given in [34], 𝐷 = 𝐴𝜎𝐵∗, where 𝜎 is (𝑛 + 1) × (𝑛) matrix given in [35] as 𝜎 = [0 0 0 … 0 1 0 0 … 0 0 2 0 … 0 ⋮ ⋮ ⋮ ⋱ 0 0 0 0 … 𝑛 ] and, 𝐵∗ is (𝑛) × (𝑛 + 1) matrix define as 𝐵∗ = [𝐴[1] −1 𝐴[2] −1 ⋮ 𝐴[𝑛] −1 ] , 𝐴[𝑘] −1 𝑖𝑠 𝑘𝑡ℎ 𝑟𝑜𝑤 𝑜𝑓 𝐴−1, 𝑘 = 1,2, ⋯ , 𝑛. We can generalize Eq. (5) as follows ⅆ𝑛 ⅆ𝑥𝑛 𝛷(𝑥) = (𝐷)𝑛𝛷(𝑥), 𝑤ℎ𝑒𝑟𝑒 𝑛 = 1, 2, . . .. (6) Let us now introduce the operational matrix of the fractional derivatives 𝐷(𝛼) , 𝛼 > 0, of size (𝑛 + 1) × (𝑛 + 1) is defined by the concept of Caputo in [35] 𝐷(𝛼) = [∑𝑛𝑗=⌈𝛼⌉   𝜔0,𝑗,0 ∑ 𝑛 𝑗=⌈𝛼⌉   𝜔0,𝑗,1 … ∑ 𝑛 𝑗=⌈𝛼⌉   𝜔0,𝑗,𝑛 ⋮ ⋮ ⋮ ∑𝑛𝑗=⌈𝛼⌉   𝜔𝑖,𝑗,0 ∑ 𝑛 𝑗=⌈𝛼⌉   𝜔𝑖,𝑗,1 … ∑ 𝑛 𝑗=⌈𝛼⌉   𝜔𝑖,𝑗,𝑛 ⋮ ⋮ ⋮ ∑𝑛𝑗=⌈𝛼⌉   𝜔𝑛,𝑗,0 ∑ 𝑛 𝑗=⌈𝛼⌉   𝜔𝑛,𝑗,1 … ∑ 𝑛 𝑗=⌈𝛼⌉   𝜔𝑛,𝑗,𝑛 ] here 𝜔𝑖,𝑗,𝑙 is given by: 𝜔𝑖,𝑗,𝑙 = (−1) 𝑗−𝑖 (𝑛 𝑖 )(𝑛 − 𝑖 𝑗 − 𝑖 ) 𝛤(𝑗 + 1) 𝛤(𝑗 + 1 − 𝑞) ∑ 𝑛 𝑘=0   𝜆𝑙𝑘 𝜇𝑘𝑗 , where 𝜆𝑙𝑘 is given in Eq. (4) and IHJPAS. 36 (3) 2023 386 𝜇 𝑘𝑗 = ∑ 𝑛 𝑠=𝑘   (−1)𝑠−𝑘 (𝑛 𝑘 )(𝑛 − 𝑘 𝑠 − 𝑘 ) 1 𝑗 − 𝛼 + 𝑠 + 1 . 3. 3 The steps of applying the BOM method for solving the Delay-Pantograph equation In this paper, we consider the fractional neutral pantograph differential equation 𝐷𝛼 𝑦(𝑥) = 𝑎(𝑥)𝑦(𝑝𝑥) + 𝑏(𝑥)𝐷𝛾 𝑦(𝑝𝑥) + 𝑑(𝑥)𝑦(𝑥) + 𝑔(𝑥), (7) where 𝑥 ∈ [0,1], 0 < 𝛾 ≤ 𝛼 ≤ 2, 0 < 𝑝 < 1. 𝑦(𝑥) is the unknown function, 𝑎, 𝑏, 𝑑, 𝑔 ∈ 𝐶[0,1], with initial and boundary conditions. Now to apply the BOM method, we follow the steps below: I- We will approximate unknown function 𝑦 (𝑥) by Bernstein polynomials as (3), and we have 𝐷𝑦(𝑥) = 𝐶𝑇 𝐷 𝛷(𝑥), 𝐷𝑦(𝑝𝑥) = 𝐶𝑇 𝐷 𝛷(𝑝𝑥),…, 𝐷𝑛 𝑦(𝑥) = 𝐶𝑇 𝐷𝑛 𝛷(𝑥), 𝐷𝑛 𝑦(𝑝𝑥)𝐶𝑇 𝐷𝑛 𝛷(𝑝𝑥), 𝐷𝛼 𝑦(𝑥) = 𝐶𝑇 𝐷(𝛼)𝛷(𝑥), 𝐷𝛼 𝑦(𝑝𝑥) = 𝐶𝑇 𝐷(𝛼)𝛷(𝑝𝑥), 𝑤ℎ𝑒𝑟𝑒 𝛼 𝑖𝑠 𝑜𝑟𝑑𝑒𝑟 𝑜𝑓 𝑡ℎ𝑒 𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛𝑎𝑙 𝑑𝑒𝑟𝑖𝑣𝑎𝑡𝑖𝑣𝑒. II- We will substitute all initial or boundary conditions given in the problem. III- We use the roots of Chebyshev polynomials to help reduce interpolation errors as the collocation node, 𝑥𝑖 = 1 2 + 1 2𝑐𝑜𝑠 ((2𝑖 + 1) 𝜋 2𝑛 ) , 𝑖 = 0,1, … , 𝑛 − 1. IV- We substitute the approximate polynomial 𝑦(𝑥) and its derivatives in Eq. (7) to get the system of algebraic equations that we can solve using computer programs like Mathematica and thus get the vector values 𝐶𝑇 . As a result, we will get an approximate solution to the problem. 4. The shifted Legendre polynomials and their operational matrix [36] Shifted Legendre polynomials will be discussed in this Section. Then we will describe the method (LOM) and how to use it to solve the problem. 4. 1. Legendre polynomials Legendre polynomials are known for the interval [-1,1] and can be calculated using the regression relationship as following 𝐿𝑖+1(𝑧) = 2𝑖 + 1 𝑖 + 1 𝑧𝐿𝑖 (𝑧) − 𝑖 𝑖 + 1 𝐿𝑖−1(𝑧), 𝑖 = 1,2, … (8) where 𝐿0(𝑧) = 1 𝑎𝑛𝑑 𝐿1(𝑧) = 𝑧. Substitute the variable 𝑧 = 2𝑥 − 1 for these polynomials in the interval [0,1], and create the so- called shifted Legendre polynomial, let 𝑃𝑖 (𝑥) be a symbol for the shifted Legendre polynomial. 𝑃𝑖 (𝑥) is calculated as follows: IHJPAS. 36 (3) 2023 387 𝑃𝑖+1(𝑥) = (2𝑖 + 1)(2𝑥 − 1) (𝑖 + 1) 𝑃𝑖 (𝑥) − 𝑖 𝑖 + 1 𝑃𝑖−1(𝑥) , 𝑖 = 1,2, … , (9) here 𝑃0(𝑥) = 1 and 𝑃1(𝑥) = 𝑥. The following formulations of the shifted Legendre polynomial of degree 𝑖 analytic form are found in the source [36]: 𝑃𝑖 (𝑥) = ∑ 𝑖 𝑘=0   (−1)𝑖+𝑘 (𝑖 + 𝑘)! (𝑖 − 𝑘)! 𝑥𝑘 (𝑘!)2 . (10) Whereas: 𝑃𝑖 (1) = 1 𝑎𝑛𝑑 𝑃𝑖 (0) = (−1) 𝑖. The requirement of orthogonality is: ∫ 1 0   𝑃𝑖 (𝑥)𝑃𝑗 (𝑥)𝑑𝑥 = { 1 2𝑖+1 𝑓𝑜𝑟 𝑖 = 𝑗, 0 𝑓𝑜𝑟 𝑖 ≠ 𝑗. The shifted Legendre polynomials' power series is written as follows: 𝑃𝑖 (𝑥) = ∑ 𝑖 𝑘=0   (−1)𝑖+𝑘 (𝑖 + 𝑘 𝑘 )(𝑖 𝑘 )𝑥𝑘 . Now we can approximate any function based on shifted Legendre polynomial as follows: 𝑦(𝑥) = ∑ ∞ 𝑗=0   𝑐𝑗 𝑃𝑗 (𝑥), where 𝑐𝑗 = (2𝑗 + 1) ∫ 1 0  𝑦(𝑥)𝑃𝑗 (𝑥)𝑑𝑥, 𝑗 = 1,2, … Only the first (𝑛 + 1) terms of shifted Legendre polynomials are considered in practice. Next, there is 𝑦(𝑥) = ∑ 𝑛 𝑗=0   𝑐𝑗 𝑃𝑗 (𝑥) = 𝐶𝑇 𝛷(𝑥), (11) where 𝐶𝑇 = [𝑐0, … , 𝑐𝑛], 𝑡ℎ𝑒 𝑠ℎ𝑖𝑓𝑡𝑒𝑑 𝐿𝑒𝑔𝑒𝑛𝑑𝑟𝑒 𝑐𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 𝑣𝑒𝑐𝑡𝑜𝑟 , 𝛷(𝑥) = [𝑃0(𝑥), 𝑃1(𝑥), … , 𝑃𝑛 (𝑥)] 𝑇 𝑡ℎ𝑒 𝑠ℎ𝑖𝑓𝑡𝑒𝑑 𝐿𝑒𝑔𝑒𝑛𝑑𝑟𝑒 𝑣𝑒𝑐𝑡𝑜𝑟. 4. 2. Operational Matrix of Fractional Derivative for Shifted Legendre Polynomials (LOM) The operational matrix of derivatives defined based on shifted Legendre polynomials is denoted by 𝐷(1) of degree (𝑛 + 1) × (𝑛 + 1) and can be expressed in the following form [36] IHJPAS. 36 (3) 2023 388 𝐷(1) = (𝑑𝑖𝑗 ) = {2(2𝑗 + 1), 𝑓𝑜𝑟 𝑗 = 𝑖 − 𝑘, 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 {𝑘 = 1,3, … , 𝑛, 𝑖𝑓 𝑛 𝑜𝑑𝑑. 𝑘 = 1,3, … , 𝑛 − 1, 𝑖𝑓 𝑛 𝑒𝑣𝑒𝑛. The vector's derivative 𝛷(𝑥) can be written as follows 𝑑𝛷(𝑥) 𝑑𝑥 = 𝐷(1)𝛷. (12) We can generalize Equation (12) as 𝑑𝑛 𝛷(𝑥) 𝑑𝑥𝑛 = (𝐷(1)) 𝑛 𝛷(𝑥), 𝑛 = 1, 2, . . . . . . (13) Let us now represent the fractional derivative operational matrix 𝐷(𝛼), 𝛼 > 0 , of size (𝑛 + 1) × (𝑛 + 1) defined by the Caputo concept in the source [36] as follows: 𝐷𝛼 = [0 0 ⋯ 0 ⋮ ⋮ ⋯ ⋮ 0 0 ⋯ 0 ∑ ⌈𝛼⌉ 𝑘=⌈𝛼⌉   𝜃⌈𝛼⌉,0,𝑘 ∑ ⌈𝛼⌉ 𝑘=⌈𝛼⌉   𝜃⌈𝛼⌉,1,𝑘 ⋯ ∑ ⌈𝛼⌉ 𝑘=⌈𝛼⌉   𝜃⌈𝛼⌉,𝑛,𝑘 ⋮ ⋮ ⋯ ⋮ ∑𝑖𝑘=⌈𝛼⌉   𝜃𝑖,0,𝑘 ∑ 𝑖 𝑘=⌈𝛼⌉   𝜃𝑖,1,𝑘 ⋯ ∑ 𝑖 𝑘=⌈𝛼⌉   𝜃𝑖,𝑛,𝑘 ⋮ ⋮ ⋯ ⋮ ∑𝑛𝑘=⌈𝛼⌉   𝜃𝑛,0,𝑘 ∑ 𝑛 𝑘=⌈𝛼⌉   𝜃𝑛,1,𝑘 ⋯ ∑ 𝑛 𝑘=⌈𝛼⌉   𝜃𝑛,𝑛,𝑘 ], where 𝜃𝑖,𝑗,𝑘 𝑖𝑠 𝑑𝑒𝑓𝑖𝑛𝑒𝑑 𝑏𝑦 𝜃𝑖,𝑗,𝑘 = (2𝑗 + 1) ∑ 𝑗 𝑙=0   (−1)𝑖+𝑗+𝑘+𝑙 (𝑖 + 𝑘)! (𝑙 + 𝑗)! (𝑖 − 𝑘)! 𝑘! 𝛤(𝑘 − 𝛼 + 1)(𝑗 − 𝑙)! (𝑙!)2(𝑘 + 𝑙 − 𝛼 + 1) . It is worth noting that the initial ⌈𝛼⌉ rows in 𝐷(𝛼)are all zero. 4.3 The steps of applying the LOM method for solving the Delay-Pantograph equation I- We will approximate unknown function 𝑦 (𝑥) by Legendre polynomials as Eq. (11), and we have 𝐷𝑦(𝑥) = 𝐶𝑇 𝐷 𝛷(𝑥), 𝐷𝑦(𝑝𝑥) = 𝐶𝑇 𝐷 𝛷(𝑝𝑥),…, 𝐷𝑛 𝑦(𝑥) = 𝐶𝑇 𝐷𝑛 𝛷(𝑥), 𝐷𝑛 𝑦(𝑝𝑥)𝐶𝑇 𝐷𝑛 𝛷(𝑝𝑥), 𝐷𝛼 𝑦(𝑥) = 𝐶𝑇 𝐷(𝛼)𝛷(𝑥), 𝐷𝛼 𝑦(𝑝𝑥) = 𝐶𝑇 𝐷(𝛼)𝛷(𝑝𝑥), 𝑤ℎ𝑒𝑟𝑒 𝛼 𝑖𝑠 𝑜𝑟𝑑𝑒𝑟 𝑜𝑓 𝑡ℎ𝑒 𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛𝑎𝑙 𝑑𝑒𝑟𝑖𝑣𝑎𝑡𝑖𝑣𝑒. II- We will substitute all initial or boundary conditions given in the problem. III- We use the roots of Chebyshev polynomials can help reduce interpolation errors as the collocation node, 𝑥𝑖 = 1 2 + 1 2𝑐𝑜𝑠 ((2𝑖 + 1) 𝜋 2𝑛 ) , 𝑖 = 0,1, … , 𝑛 − 1. IV- We substitute the approximate polynomial 𝑦(𝑥) and its derivatives in Eq. (7) to get the system of algebraic equations that we can solve using computer programs like Mathematica and thus get the vector values 𝐶𝑇 as a result, we will get an approximate solution to the problem. 5. Illustrative examples IHJPAS. 36 (3) 2023 389 In this Section, some numerical examples are given to demonstrate the applicability and accuracy of the proposed method. 𝑀𝑎𝑡ℎ𝑒𝑚𝑎𝑡𝑖𝑐𝑎® 12 is used for all numerical calculations. Example 1: - Consider the following fractional neutral pantograph differential equation 𝐷𝛼 𝑦(𝑥) = −𝑦(𝑥) + 0.1𝑦(0.8𝑥) + 0.5𝐷𝛼 𝑦(0.8𝑥) + 𝑔(𝑥), 𝑥 ∈ (0,1], (14) 𝑦(0) = 0, 𝑤𝑖𝑡ℎ 0 < 𝛼 ≤ 1. Now to solve Eq. (14) by (BOM), the exact solution to this problem for 𝛼 = 0.8 is 𝑦(𝑥) = 𝑥3.8 with (𝑥) = 2.211894885744887𝑥3 + 0.9571706039258614𝑥3.8 , and with 𝑛 = 2, by applying the suggested method in Section (3-2), we have 𝐷0.8 = [−1.14334 − 1.55 − 0.276996 1.17281 0.872243 − 1.55 − 0.0294677 0.677756 1.82699 ] By applying the condition given in the problem and following the solution steps mentioned in Section )3-3), a system of algebraic equations will be produced, and by solving this system, we get the values of 𝑐0 = 0, 𝑐1 = 0.0118614, 𝑐2 = 0.472083, the approximate solution is 𝑦(𝑥) = 0.0237227𝑥 + 0.44836𝑥2. If 𝑛 = 15, we have the following approximate solution: 𝑦(𝑥) = 9.543602968653811 × 10−17 + 0.000006511676397407056𝑥 − 0.0007034719961872769𝑥2 + 0.03467000244609836𝑥3 + 1.4120584194534938𝑥4 − 1.2900408484650545𝑥5 + 2.1428840263535394𝑥6 − 2.118319283516632𝑥7 − 0.6752474254279832𝑥8 + 6.081788398525532𝑥9 − 10.489702105988727𝑥10 + 10.517845113990916𝑥11 − 7.003103640066911𝑥12 + 3.1968135804363556𝑥13 − 0.9498166132520964𝑥14 + 0.14086738703346668𝑥15. Table 1 shows the absolute error|𝑦𝑒𝑥𝑎𝑐𝑡 − 𝑦𝑎𝑝𝑝.| value for the some 𝑥 ∈ (0,1] at 𝑛= 9, 11, 13, and 15. Table 1: The absolute error of Ex. 1 for 𝑛=9,11,13,15 by (BOM) 𝒙 𝒏 = 𝟗 𝒏 = 𝟏𝟏 𝒏 = 𝟏𝟑 𝒏 = 𝟏𝟓 0.1 1.21756×10-7 2.92165×10-8 6.77509×10-8 3.19953×10-8 0.2 1.07484×10-7 3.66589×10-8 7.85818×10-8 3.32075×10-8 0.3 1.40511×10-7 5.90616×10-8 5.80438×10-8 2.87904×10-8 0.4 3.05621×10-7 6.14799×10-8 2.5095×10-8 1.82177×10-8 0.5 5.54154×10-8 2.54537×10-8 9.48675×10-9 4.38442×10-9 0.6 2.92554×10-7 3.27859×10-8 2.82086×10-8 6.08513×10-10 0.7 3.1252×10-7 8.2099×10-8 3.73298×10-8 1.17502×10-8 0.8 6.18788×10-8 1.2187×10-7 3.22221×10-8 1.2558×10-9 0.9 2.07441×10-7 1.90066×10-7 1.81549×10-8 1.60439×10-8 In order to use the LOM approach to solve example 1, we must follow the steps outlined in section (4-3). If we assume 𝑛 = 2, we get 𝐷0.8 = [0 0 0 1.81521 0.495057 − 0.206274 − 0.495057 4.08422 1.06084 ] The approximate solution 𝑦(𝑥) = 0.0237227𝑥 + 0.44836𝑥2, and the value of 𝑐0 = 0.161315, 𝑐1 = 0.236041, 𝑐2 = 0.0747266. If 𝑛 = 15, the approximate solution 𝑦(𝑥) = −5.917880389286895 × 10−17 + 0.000004941255572793077𝑥 − 0.0005604595883997804𝑥2 + 0.030640435190110164𝑥3 + 1.4662615402598584𝑥4 − 1.7082026330161135𝑥5 + 4.183434690591156𝑥6 − 8.78814674633918𝑥7 + IHJPAS. 36 (3) 2023 390 14.437646065525456𝑥8 − 18.130978138638397𝑥9 + 17.216360479207815𝑥10 − 12.200969372917305𝑥11 + 6.3079153582312015𝑥12 − 2.277429378098255𝑥13 + 0.5224877022689723𝑥14 − 0.058464614185574325𝑥15. We present Table 2, in which we will show the absolute error value of 𝑛=9,11,13, and 15 Table 2: The absolute error of Ex. 1 for 𝑛=9,11,13,15 by (LOM) 𝒙 𝒏 = 𝟗 𝒏 = 𝟏𝟏 𝒏 = 𝟏𝟑 𝒏 = 𝟏𝟓 0.1 1.21719×10-7 2.88811×10-8 2.34826×10-9 1.72445×10-8 0.2 1.07453×10-7 3.60137×10-8 1.13909×10-8 1.44977×10-8 0.3 1.40539×10-7 5.78547×10-8 9.88362×10-9 1.44738×10-8 0.4 3.05492×10-7 6.08407×10-8 9.05138×10-9 1.08162×10-8 0.5 5.55965×10-8 2.58773×10-8 3.71809×10-8 1.03386×10-9 0.6 2.92938×10-7 3.09694×10-8 6.09679×10-8 8.16747×10-9 0.7 3.1172×10-7 8.21428×10-8 6.37435×10-8 1.17671×10-8 0.8 6.3123×10-8 1.21236×10-7 5.67603×10-8 5.21886×10-9 0.9 2.09417×10-7 1.93688×10-7 6.72274×10-8 5.09303×10-9 In Figure 1, we show a plot of both the exact solution and the approximate solution of the two proposed methods for 𝑛=15. Figure 1: The comparison of the exact solution and the approximate solution of (BOM, LOM) methods for Ex. 1 at 𝑛=15 Also, in Figure 2, logarithmic plots of the greatest absolute error of the approximate solution are given using Bernstein or Legendre polynomials with operational matrices to solve example 1 from 𝑛 = 2 to 𝑛 = 15. IHJPAS. 36 (3) 2023 391 Figure 2: Logarithmic plots for the 𝑀𝐴𝐸𝑛 of Ex. 1 form 𝑛=2 to 𝑛=15 We note from the above results that the two proposed methods are efficient because their results are highly accurate compared to the results of the exact solution, as shown in Figure 1. Also, we can see the efficiency of the two proposed methods in Figure 2, which shows the decrease in absolute error when the value of 𝑛 is increased. In order to show which of the two methods is more efficient, we present Table 3, in which we will show the absolute error values of the two proposed methods at 𝑛 = 15 Table 3: The absolute error |𝑦𝑒𝑥𝑎𝑐𝑡 − 𝑦𝑎𝑝𝑝.| of (BOM, LOM) for Ex. 1 at 𝑛=15 𝒙 Absolute error (B) Absolute error (L) 0.1 3.19953×10-8 1.72445×10-8 0.2 3.32075×10-8 1.44977×10-8 0.3 2.87904×10-8 1.44738×10-8 0.4 1.82177×10-8 1.08162×10-8 0.5 4.38442×10-9 1.03386×10-9 0.6 6.08513×10-10 8.16747×10-9 0.7 1.17502×10-8 1.17671×10-8 0.8 1.2558×10-9 5.21886×10-9 0.9 1.60439×10-8 5.09303×10-9 We notice from Table 3 that the values of Legendre's absolute error are less than those of Bernstein’s absolute error, which indicates that (LOM) method is more efficient than (BOM) method in solving example 1. Example 2: Consider the following fractional pantograph equation 𝐷2𝑦(𝑥) + 𝐷 3 2𝑦(𝑥) + 𝑦(𝑥) = 𝑦 ( 𝑥 2 ) + 3𝑥2 4 + 4√ 𝑥 𝜋 + 2, 𝑥 ∈ [0,1], (15) subject to 𝑦(0) = 0, 𝑦(1) = 1, and the exact solution to this problem is 𝑦(𝑥) = 𝑥2. When applying the (BOM) method for solving Eq. (15) at 𝑛 = 2, we have 𝐷 3 2 = [ 24 35√𝜋 24 7√𝜋 136 35√𝜋 − 48 35√𝜋 − 48 7√𝜋 − 272 35√𝜋 24 35√𝜋 24 7√𝜋 136 35√𝜋 ], We also obtained the values of the unknown vector 𝐶𝑇 , which are as follows 𝑐0 = 1.18599 × 10−17, 𝑐1 = 0.00228219, 𝑐2 = 1. IHJPAS. 36 (3) 2023 392 We also come up with the following approximate solution: 𝑦(𝑥) = 1.18599 × 10−17 + 0.00456438𝑥 + 0.995436𝑥2. If 𝑛 = 20, we have the following approximate solution: 𝑦(𝑥) = −2.786853970684164 × 10−15 − 0.000030030531854295995𝑥 + 1.006670875743893𝑥2 − 0.31622314969631304𝑥3 + 7.39417733179431𝑥4 − 101.22636028238139𝑥5 + 882.5852080762097𝑥6 − 5099.047105444493𝑥7 + 19465.37892066044𝑥8 − 44677.52886989678𝑥9 + 27582.502303316025𝑥10 + 211915.94892108513𝑥11 − 939140.5398020139𝑥12 + 2196019.891017601𝑥13 − 3463946.68783369𝑥14 + 3900531.650722769𝑥15 − 3162032.689661729𝑥16 + 1810684.418812781𝑥17 − 696899.7151342098𝑥18 + 162020.35868317497𝑥19 − 17212.384417226098𝑥20. Table 4 shows the absolute error|𝑦𝑒𝑥𝑎𝑐𝑡 − 𝑦𝑎𝑝𝑝.| value for the some 𝑥 ∈ [0,1] at 𝑛= 14,16,18, and 20 Table 4: The absolute error of Ex. 2 for 𝑛=14,16,18,20 by (BOM) 𝒙 𝒏 = 𝟏𝟒 𝒏 = 𝟏𝟔 𝒏 = 𝟏𝟖 𝒏 = 𝟐𝟎 0.1 3.3016×10-6 2.0304×10-6 1.679×10-6 1.42028×10-6 0.2 5.64738×10-6 3.96901×10-6 2.44675×10-6 1.92033×10-6 0.3 5.90691×10-6 4.01465×10-6 3.34437×10-6 2.143×10-6 0.4 6.08868×10-6 4.68942×10-6 2.963×10-6 2.48202×10-6 0.5 6.6422×10-6 3.86708×10-6 3.25554×10-6 2.06817×10-6 0.6 4.61383×10-6 4.11001×10-6 2.54331×10-6 1.97996×10-6 0.7 4.97604×10-6 2.68729×10-6 2.29425×10-6 1.58232×10-6 0.8 2.62815×10-6 2.43508×10-6 1.62869×10-6 9.01401×10-7 0.9 2.04652×10-6 1.1369×10-6 6.75106×10-7 3.45435×10-7 To solve example 2 using the (LOM) technique, we must follow the steps indicated in section (4- 3). If we take n=2, we get 𝐷 3 2 = [0 0 0 0 0 0 16 √𝜋 48 5√𝜋 − 16 7√𝜋 ], We also obtain the values of the unknown vector 𝐶𝑇 , which are as follows 𝑐0 = 0.334094, 𝑐1 = 0.5, 𝑐2 = 0.165906. and their approximate solution is 𝑦(𝑥) = 2.77555756 × 10−17 + 0.00456438𝑥 + 0.995436𝑥2. If 𝑛 = 20, the approximate solution will be as follows: 𝑦(𝑥) = −7.827619401012162 × 10−17 − 0.000030068562737359008𝑥 + 1.0066701273460272𝑥2 − 0.3161559593254199𝑥3 + 7.391520839002448𝑥4 − 101.1675850252221𝑥5 + 881.7639994208571𝑥6 − 5091.2533573008905𝑥7 + 19412.68134451585𝑥8 − 44415.30443366349𝑥9 + 26600.722236690082𝑥10 + 214722.2039331961𝑥11 − 945315.3904772103𝑥12 + 2206507.388919022𝑥13 − 3477652.419776158𝑥14 + 3914177.755634204𝑥15 − 3172189.258528861𝑥16 + 1816152.725780014𝑥17 − 698910.3552130745𝑥18 + 162471.92620954153𝑥19 − 17259.100690250787𝑥20. IHJPAS. 36 (3) 2023 393 Table 5 shows the absolute error|𝑦𝑒𝑥𝑎𝑐𝑡 − 𝑦𝑎𝑝𝑝.| value for the some 𝑥 ∈ [0,1] at 𝑛= 14,16,18, and 20 Table 5: The absolute error of Ex. 2 for 𝑛=14,16,18,20 by (LOM) 𝒙 𝒏 = 𝟏𝟒 𝒏 = 𝟏𝟔 𝒏 = 𝟏𝟖 𝒏 = 𝟐𝟎 0.1 3.30159×10-6 2.03041×10-6 1.67929×10-6 1.41633×10-6 0.2 5.64737×10-6 3.96904×10-6 2.44733×10-6 1.91244×10-6 0.3 5.9069×10-6 4.0147×10-6 3.34525×10-6 2.13118×10-6 0.4 6.08866×10-6 4.68948×10-6 2.96416×10-6 2.4665×10-6 0.5 6.64218×10-6 3.86716×10-6 3.25698×10-6 2.04886×10-6 0.6 4.6138×10-6 4.11011×10-6 2.54504×10-6 1.95697×10-6 0.7 4.97601×10-6 2.6874×10-6 2.29617×10-6 1.55857×10-6 0.8 2.62812×10-6 2.43526×10-6 1.63045×10-6 8.94486×10-7 0.9 2.04649×10-6 1.13714×10-6 6.78982×10-7 4.05584×10-7 In Figure 3, we show a plot of both the exact solution and the approximate solution of the two proposed methods of example 2 for 𝑛=20. Figure 3: The comparison of the exact solution and the approximate solution of (BOM, LOM) methods for Ex. 2 at 𝑛=20 In Figure 4, logarithmic plots of the greatest absolute error of the approximate solution are given using Bernstein or Legendre polynomials with operational matrices to solve example 2 from 𝑛 = 2 to 𝑛 = 20. IHJPAS. 36 (3) 2023 394 Figure 4: Logarithmic plots for the 𝑀𝐴𝐸𝑛 of Ex. 2 form 𝑛=2 to 𝑛=20 The two proposed approaches are efficient, as seen by the above results because their results are accurate when compared to the precise solution's results, as shown in Figure 1. Figure 2 demonstrates the decrease in absolute error when the value of 𝑛 is increased, demonstrating the efficacy of the two proposed techniques. Table 6 illustrates the absolute error numbers for (BOM, LOM) methods at 𝑛 = 20 to determine which of the two strategies is superior for solving example 2. Table 6: The absolute error |𝑦𝑒𝑥𝑎𝑐𝑡 − 𝑦𝑎𝑝𝑝.| of (BOM, LOM) for Ex. 2 at 𝑛=20 𝒙 Absolute error (B) Absolute error (L) 0.1 1.42028×10 -6 1.41633×10-6 0.2 1.92033×10 -6 1.91244×10-6 0.3 2.143×10 -6 2.13118×10-6 0.4 2.48202×10 -6 2.4665×10-6 0.5 2.06817×10 -6 2.04886×10-6 0.6 1.97996×10 -6 1.95697×10-6 0.7 1.58232×10 -6 1.55857×10-6 0.8 9.01401×10 -7 8.94486×10-7 0.9 3.45435×10 -7 4.05584×10-7 We notice from the above table that the greatest absolute error we acquired in the method (BOM) is 2.48202×10-6, but in the method (LOM) is 2.4665×10-6, indicating that the method (LOM) is the best in solving Ex. 2. 6. Conclusion In this study, we propose the methods of operational matrices that depend on Bernstein polynomials and Shifted Legendre polynomials for fractional derivatives in solving fractional delay pantograph equations. The numerical results of the approximate solutions in examples 1 and 2 were given accuracy when compared with the exact solutions, proving the effectiveness and efficiency of the proposed methods. Figures 2 and 4 showed the absolute error of the two examples for more than 𝑛, where when the values of 𝑛 increase, the value of the absolute error decreases, which shows the accuracy and success of the proposed methods. Tables 3 and 6 show that the method (LOM) is the best solution for the two examples. Mathematica®12 program was used to find the numerical results. IHJPAS. 36 (3) 2023 395 References 1. Jamil. B.; Anwar. M.S.; Rasheed. A.; Irfan. M. MHD Maxwell flow modeled by fractional derivatives with chemical reaction and thermal radiation. Chinese Journal of Physics. 2020, 67, 512-533, doi: org/10.1016/j.cjph.2020.08.012. 2. Alzubaidi. W.K.; Shaker. S.H. Methods of Secure Routing Protocol in Wireless Sensor Networks. Journal of AL-Qadisiyah for computer science and mathematics. 2018, 10, 3, 3855- , doi: 10.29304/jqcm.2018.10.3.437. 3. Alchikh. R.; Khuri. S.A. Numerical solution of a fractional differential equation arising in optics. Optik. 2020, 208, 163911 –163919, doi: org/10.1016/j.ijleo.2019.163911. 4. Bahaa. G. M. Optimal control problem for variable-order fractional differential systems with time delay involving Atangana–Baleanu derivatives. Chaos, Solitons, Fractals. 2019, 122, 129- 142, doi: org/10.1016/j.chaos.2019.03.001. 5. Tarasov. V.E.; Aifantis. E. C. On fractional and fractal formulations of gradient linear and nonlinear elasticity. Acta Mechanica. 2019, 230, 6, 2043-2070, doi: org/10.1007/s00707-019- 2373-x. 6. Goufo. E.F.D.; Nieto. J.J. Attractors for fractional differential problems of transition to turbulent flows. Journal of Computational and Applied Mathematics. 2018, 339, 329-342, doi: org/10.1016/j.cam.2017.08.026. 7. Mechee. M. S.; Al-Shaher. O. I.; Al-Juaifri. G. A. Haar wavelet technique for solving fractional differential equations with an application. In AIP Conference Proceedings. 2019, 2086(1), 030025, doi:org/10.1063/1.5095110. 8. Agarwal. P.; El-Sayed. A. A. Non-standard finite difference and Chebyshev collocation methods for solving fractional diffusion equation. Physica A: Statistical Mechanics and Its Applications. 2018, 500, 40-49, doi: https://doi.org/10.1016/j.physa.2018.02.014. 9. Alizadeh. A.; Effati, S. Modified Adomian decomposition method for solving fractional optimal control problems. Transactions of the Institute of Measurement and Control. 2018, 40(6), 2054- 206, doi: https://doi.org/10.1177/0142331217700243. 10. Thabet. H.; Kendre. S. New modification of Adomian decomposition method for solving a system of nonlinear fractional partial differential equations. Int. J. Adv. Appl. Math. Mech. 2019, 6(3), 1-13, https://www.researchgate.net/publication/332014078. 11. Ali. K. A., Approximate solution for fuzzy differential algebraic equations of fractional order using Adomain decomposition method. Ibn Al-Haitham J. for Pure&Appl. Sci.2017, 30(2), 202-213. 12. Javeed. S.; Baleanu. D.; Waheed. A.; Shaukat Khan. M.; Affan. H. Analysis of homotopy perturbation method for solving fractional order differential equations. Mathematics. 2019, 7(1), 40, doi: https://doi.org/10.3390/math7010040. 13. Hamoud. A.; Ghadle. K. Usage of the homotopy analysis method for solving fractional Volterra-Fredholm integro-differential equation of the second kind. Tamkang Journal of Mathematics. 2018, 49(4), 301-315, doi: https://doi.org/10.5556/j.tkjm.49.2018.2718. 14. Jiang. Y.; Xu. X. A monotone finite volume method for time fractional Fokker-Planck equations. Science China Mathematics. 2019, 62(4), 783-794, doi: https://doi.org/10.1007/s11425-017-9179-x . 15. Baleanu. D.; Jassim. H. K.; Khan. H. A modification fractional variational iteration method for solving non-Linear gas dynamic and coupled Kdv equations involving local fractional operators. 2018, http://hdl.handle.net/20.500.12416/2276. about:blank about:blank about:blank about:blank https://doi.org/10.1016/j.ijleo.2019.163911 about:blank about:blank about:blank about:blank about:blank about:blank https://doi.org/10.1016/j.physa.2018.02.014 https://doi.org/10.1177/0142331217700243 https://www.researchgate.net/publication/332014078 https://doi.org/10.3390/math7010040 https://doi.org/10.5556/j.tkjm.49.2018.2718 https://doi.org/10.1007/s11425-017-9179-x http://hdl.handle.net/20.500.12416/2276 IHJPAS. 36 (3) 2023 396 16. Raslan. K. R.; Ali. K. K.; Mohamed. E. M. Spectral Tau method for solving general fractional order differential equations with linear functional argument. Journal of the Egyptian Mathematical Society. 2019, 27(1), 1-16, doi: https://doi.org/10.1186/s42787-019-0039-4. 17. Kuang. Y. (Ed.); Delay differential equations: with applications in population dynamics. Academic Press. 1993. 18. MacDonald. N.; MacDonald. N. Biological delay systems: linear stability theory. Cambridge University Press. 2008, www.cambridge.org/9780521340847. 19. Dehghan. M.; Shakeri, F. The use of the decomposition procedure of Adomian for solving a delay differential equation arising in electrodynamics. Physica Scripta. 2008, 78(6), 065004, doi:https://doi.org/10.1088/0031-8949/78/06/065004. 20. Magin. R.L. Fractional calculus models of complex dynamics in biological tissues. Comput Math Appl. 2010, 59,1586–1593, doi: https://doi.org/10.1016/j.camwa.2009.08.039. 21. Aiello. W. G.; Freedman. H. I.; Wu. J. Analysis of a model representing stage-structured population growth with state-dependent time delay. SIAM Journal on Applied Mathematics. 1992, 52(3), 855-869, doi: https://doi.org/10.1137/0152048 22. Rabiei. K.; Ordokhani. Y. Solving fractional pantograph delay differential equations via fractional-order Boubaker polynomials. Engineering with Computers. 2019, 35(4), 1431- 1441. doi:https://doi.org/10.1007/s00366-018-0673-8. 23. Yuttanan. B.; Razzaghi. M.; Vo. T. N. A fractional‐order generalized Taylor wavelet method for nonlinear fractional delay and nonlinear fractional pantograph differential equations. Mathematical Methods in the Applied Sciences.2021, 44(5), 4156-4175, doi: https://doi.org/10.1002/mma.7020. 24. Avazzadeh. Z.; Heydari. M. H.; Mahmoudi. M. R. An approximate approach for the generalized variable-order fractional pantograph equation. Alexandria Engineering Journal. 2020, 59(4), 2347-2354., doi: https://doi.org/10.1016/j.aej.2020.02.028. 25. Sedaghat. S.; Ordokhani. Y.; Dehghan. M. Numerical solution of the delay differential equations of pantograph type via Chebyshev polynomials. Communications in Nonlinear Science and Numerical Simulation. 2012, 17(12), 4815-4830, https://doi.org/10.1016/j.cnsns.2012.05.009 26. Hashemi. M. S.; Atangana. A.; Hajikhah. S. Solving fractional pantograph delay equations by an effective computational method. Mathematics and Computers in Simulation. 2020, 177, 295-305, doi: https://doi.org/10.1016/j.matcom.2020.04.026. 27. Adel. W.; Sabir, Z. Solving a new design of nonlinear second-order Lane–Emden pantograph delay differential model via Bernoulli collocation method. The European Physical Journal Plus. 2020, 135(5), 1-12, doi: https://doi.org/10.1140/epjp/s13360-020-00449-x 28. Rostamy. D.; Qasemi. S. Discontinuous Petrov-Galerkin and Bernstein-Legendre polynomials method for solving fractional damped heat-and wave-like equations. Journal of Computational and Theoretical Transport. 2015, 44(1), 1-23, doi: https://doi.org/10.1080/23324309.2014.955207. 29. Benattia. M. E.; Belghaba. K. Numerical solution for solving fractional differential equations using shifted Chebyshev wavelet. Gen. Lett. Math. 2017, 3(2), 101-110. 30. Doha. E. H.; Bhrawy. A. H.; Ezz-Eldien. S. S. A new Jacobi operational matrix: an application for solving fractional differential equations. Applied Mathematical Modelling. 2012, 36(10), 4931-4943, doi: https://doi.org/10.1016/j.apm.2011.12.031. https://doi.org/10.1186/s42787-019-0039-4 about:blank about:blank https://doi.org/10.1088/0031-8949/78/06/065004 https://doi.org/10.1016/j.camwa.2009.08.039 https://doi.org/10.1137/0152048 https://doi.org/10.1007/s00366-018-0673-8 https://doi.org/10.1002/mma.7020 https://doi.org/10.1016/j.aej.2020.02.028 https://doi.org/10.1016/j.cnsns.2012.05.009 https://doi.org/10.1016/j.matcom.2020.04.026 https://doi.org/10.1140/epjp/s13360-020-00449-x https://doi.org/10.1080/23324309.2014.955207 https://doi.org/10.1016/j.apm.2011.12.031. IHJPAS. 36 (3) 2023 397 31. Mohammadi. F.; Hosseini. M. M. A new Legendre wavelet operational matrix of derivative and its applications in solving the singular ordinary differential equations. Journal of the Franklin Institute. 2011, 348(8), 1787-1796., doi: https://doi.org/10.1016/j.jfranklin.2011.04.017 32. Ibraheem. G.H.; Al-Jawary. M.A. The operational matrix of Legendre polynomials for solving nonlinear thin film flow problems. Alexandria Engineering Journal. 2020, 59,5, 4027-4033, doi: org/10.1016/j.aej.2020.07.008. 33. Alshbool. M. H. T.; Bataineh. A. S.; Hashim. I.; Isik. O. R. Solution of fractional-order differential equations based on the operational matrices of new fractional Bernstein functions. Journal of King Saud University-Science. 2017, 29(1), 1-18, doi: https://doi.org/10.1016/j.jksus.2015.11.004. 34. Al-Jawary. M.A.; Ibraheem. G.H. Tow meshless methods for solving nonlinear ordinary differential equations in engineering and applied sciences. Nonlinear Eng. J. 2020, 9,1, 244- 255, doi: org/10.1515/nleng-2020-0012. 35. Saadatmandi. A.; Bernstein operational matrix of fractional derivatives and its applications. Applied Mathematical Modelling. 2014, 38(4), 1365-1372, doi: https://doi.org/10.1016/j.apm.2013.08.007. 36. Saadatmandi. A.; Dehghan. M. A new operational matrix for solving fractional-order differential equations. Computers Mathematics with Applications. 2010, 59(3), 1326-1336, doi:https://doi.org/10.1016/j.camwa.2009.07.006 https://doi.org/10.1016/j.jfranklin.2011.04.017 about:blank https://doi.org/10.1016/j.jksus.2015.11.004 https://doi.org/10.1515/nleng-2020-0012 https://doi.org/10.1016/j.apm.2013.08.007 https://doi.org/10.1016/j.camwa.2009.07.006