IHJPAS. 36 (4) 2023 

338 

This work is licensed under a Creative Commons Attribution 4.0 International License 

 

   

 

 
 

 

 

 

 

 

 

 

 

*Corresponding Author: othman.m.s@ihcoedu.uobaghdad.edu.iq 

 

 
  

 

 

Abctract 

This paper investigates an Effective Computational Method (ECM) based on the standard 

polynomials used to solve some nonlinear initial and boundary value problems in engineering 

and applied sciences. Moreover, the effective computational methods in this paper were improved 

by suitable orthogonal base functions, especially the Chebyshev, Bernoulli, and Laguerre 

polynomials, to obtain novel approximate solutions for some nonlinear problems. These base 

functions enable the nonlinear problem to be effectively converted into a nonlinear algebraic 

system of equations, which are then solved using Mathematica®12. The Improved Effective 

Computational Methods (I-ECMs) have been implemented to solve three applications involving 

nonlinear initial and boundary value problems: the Darcy-Brinkman-Forchheimer equation, the 

Blasius equation, and the Falkner-Skan equation, and a comparison between the proposed 

methods has been presented. Furthermore, the Maximum Error Remainder (𝑀𝐸𝑅𝑛) has been 

computed to prove the proposed methods' accuracy. The results convincingly prove that ECM 

and I-ECMs are effective and accurate in obtaining novel approximate solutions to the problems. 

 

Keywords: Darcy-Brinkman-Forchheimer equation; Blasius equation; Falkner-Skan equation; 

Chebyshev polynomials; Bernoulli polynomials; Laguerre polynomials. 

 

doi.org/10.30526/36.4.3265 

Article history: Received 03 Februray 2022, Accepted 14 March 2023, Published in October 2023 

 

 

Ibn Al-Haitham Journal for Pure and Applied Sciences 

Journal homepage: jih.uobaghdad.edu.iq 

 

Novel Approximate Solutions for Nonlinear Initial and Boundary Value Problems 

       Othman Mahdi Salih*  

Department of Mathematics, College of Education 

for Pure Sciences Ibn AL-Haitham, University of 

Baghdad, Baghdad, Iraq. 

 

 

    Majeed A. AL-Jawary  

Department of Mathematics, College of Education 

for Pure Sciences Ibn AL-Haitham, University of 

Baghdad, Baghdad, Iraq. 

 

    Mustafa Turkyilmazoglu  

-Department of Mathematics, Hacettepe 

University, Ankara, Turkey. 

-Department of Medical Research, China Medical 

University Hospital, China Medical University, 

Taichung, Taiwan. 

turkyilm@hacettepe.edu.tr 

 

 

https://creativecommons.org/licenses/by/4.0/
mailto:othman.m.s@ihcoedu.uobaghdad.edu.iq
mailto:othman.m.s@ihcoedu.uobaghdad.edu.iq
mailto:majeed.a.w@ihcoedu.uobaghdad.edu.iq
mailto:turkyilm@hacettepe.edu.tr
mailto:turkyilm@hacettepe.edu.tr


IHJPAS. 36 (4) 2023 

339 
 

1. Introduction 

There are many problems in engineering and applied sciences, such as fluid flow models, 

mechanical engineering, and mathematical physics, which can be described by nonlinear ordinary 

differential equations [1]. This leads to significant computational difficulties for nonlinear 

boundary conditions, particularly for nonlinear initial value problems [2]. Since the exact 

solutions to these problems are often complicated or occasionally may not be available. 

Therefore, there is a great need to develop efficient, novel approximate, and numerical methods 

to solve these problems [3 and 4]. 

Numerous analytical and approximate methods for solving nonlinear differential equations 

have been introduced and developed by authors around the world, such as the advanced Adomian 

decomposition method [5], the Variational Iteration Method (VIM), the Differential 

Transformation Method (DTM) [6], the finite difference methods [7], the optimal quartic B-

spline collocation method [8], the homotopy analysis method with Padé approximations [9], the 

Chebyshev operational matrix method [10], the Bernoulli matrix method [11], the Laguerre 

collocation method [12]. In particular, AL-Jawary et al. [13] have applied the Daftardar-Jafari 

Method (DJM), the Temimi-Ansari Method (TAM), and the Banach Contraction Method (BCM) 

to obtain the solution for the Jeffery-Hamel flow problem. Agom et al. [14] have implemented 

the Homotopy Perturbation Method (HPM) and the Adomian Decomposition Method (ADM) for 

solving the 12𝑡ℎ-order boundary value problems in finite domains. Also, Singh [15] used the 

modified homotopy perturbation approach to solve a set of nonlinear Lane-Emden equations. 

Ibraheem et al. [16] have recently implemented the operational matrix of Legendre 

polynomials to solve nonlinear thin-film flow problems. Gürbüz et al. [17] used the matrix 

relations between the Laguerre polynomials and their derivatives to study second-order nonlinear 

ordinary differential equations with quadratic and cubic terms and several other approximation 

methods for instance, see [18-23]. 

In recent years, approximation methods for analyzing linear systems of ordinary differential 

equations using orthogonal series have been widely developed. These are known as spectral 

methods, assuming that a truncated orthogonal series expansion can reasonably approximate the 

solution. Depending on the nature of the problem, a variety of orthogonal series have been used, 

such as the Walsh series, block-pulse, Laguerre, Chebyshev, Fourier series, and others [24]. 

Furthermore, orthogonal functions and polynomial series have attracted significant attention 

because they have been instrumental in treating various dynamical system problems. The main 

feature of this technique is that it reduces these problems to the solution of a system of algebraic 

equations by using the method of operational matrices based on orthogonal polynomials [25], 

such as Chebyshev polynomials [26], Bernoulli polynomials [27], and Laguerre polynomials 

[28], which significantly simplifies the problems and allows them to be solved by any 

computational program. 

More recently, Turkyilmazoglu [29] has proposed and used an analytical approximation 

method, namely the ECM, to solve various types of problems, such as nonlinear Lane-Emden-

Fowler equations [29], Fredholm integro-differential equations [30], Volterra-Fredholm-

Hammerstein integro-differential equations [31], heat transfer of fin problems [32], and initial 

and boundary value problems with difficult exact solutions [33]. Moreover, the approach depends 

on appropriate base functions, such as the standard polynomials. In addition, the solution of the 

nonlinear equations is transformed into a nonlinear algebraic system with unknown standard 

polynomial coefficients, which can be solved numerically or analytically with modern software. 



IHJPAS. 36 (4) 2023 

340 
 

The current paper aims to use the ECM based on the standard polynomials to solve three 

applications involving nonlinear initial and boundary value problems: the Darcy-Brinkman-

Forchheimer equation, the Blasius equation, and the Falkner-Skan equation, which are found in 

engineering and applied sciences. The main goals are to develop the ECM by introducing various 

orthogonal polynomials, such as Chebyshev, Bernoulli, and Laguerre polynomials, and to form 

a novel collection of the I-ECMs. The final goal is to implement the I-ECMs to solve these 

problems. 

The paper is organized as follows: Section two presents the mathematical formulations of three 

nonlinear models. Section three introduces the basic concepts of the proposed methods. Section 

four displays the implementation of the ECM and I-ECMs to solve three nonlinear problems and 

discusses the results. Finally, section five presents the conclusions. 

2. The Mathematical Formulations of Nonlinear Models 

2.1 The Darcy-Brinkman-Forchheimer Equation 

Consider the following steady-state, pressure-driven, fully developed parallel flow over a 

horizontal channel filled with a porous medium [34], as demonstrated in Figure 1: 

 

 

 

 

 

 

 

Figure 1. Parallel flow in a fluid-saturated porous channel [35]. 

 

The positions of the bottom and top plates are 𝑦 = ℎ and 𝑦 = − ℎ, respectively. The velocity 

takes the form 𝑢 = (𝑦(𝑥), 0, 0) and the flow is in the 𝑥-axis direction. The Darcy-Brinkman-

Forchheimer equation, which has the following form [36], is known to determine the flow in the 

channel. 

𝑦′′(𝑥) − 𝑠2 𝑦(𝑥) − 𝐹 𝑠 𝑦2(𝑥) +
1

𝑀
= 0,                                                                               (1) 

subjected to the boundary conditions: 

𝑦′(0) = 0, 𝑦(1) = 0.                                                                                                         (2) 

where 𝐹 stands for the Forchheimer number, 𝑠 for the shape parameter of the porous medium, 

and 𝑀 for the viscosity ratio. 

The Darcy-Brinkman-Forchheimer equation has been solved analytically and approximately 

using a variety of methods, such as the homotopy analysis method [37], the finite difference 

method [38], the optimal asymptotic Galerkin homotopy method [36], and the Tau homotopy 

analysis method [34]. In particular, Adewumi et al. [39] obtained the approximate solutions for 

the model by using the hybrid method in combination with the Chebyshev collocation method 

with Laplace and differential transform methods. Motsa et al. [35] implemented the spectral 

homotopy analysis approach to obtain an accurate result for the model. In addition, Abbasbandy 

et al. [40] obtained a closed-form solution of forced convection in a porous saturated channel. 

 



IHJPAS. 36 (4) 2023 

341 
 

2.2 The Blasius Equation 

The Blasius equation is the well-known third-order nonlinear ordinary differential equation 

that appeared in several boundary layer problems involving a fluid's two-dimensional laminar 

viscous flow through a flat plate. The following equation presents it as a governing equation for 

fluid dynamics [41]: 

𝑦′′′(𝑥) +
1

2
𝑦(𝑥) 𝑦′′(𝑥) = 0,                                                                                                    (3) 

subjected to the boundary conditions: 

𝑦(0) = 0,   𝑦′(0) = 0,   𝑦′(∞)  = 1.                                                                                     (4) 

The second derivative of 𝑦(𝑥) at zero is important in the Blasius equation to evaluate the shear 

stress on the plate. Numerous authors have tried to solve this problem and obtained various 

numbers for this value. For more details, see [42-44]. Therefore, the boundary conditions of the 

Blasius equation become: 

𝑦(0) = 0,   𝑦′(0) = 0,   𝑦′′(0) = 𝛼.                                                                                    (5) 

The value of 𝛼 = 0.3320573 will be utilized in the present work, as stated in [43]. 

Several numerical and analytical techniques have been used to solve the Blasius equation, such 

as the homotopy analysis method [45], the optimal homotopy asymptotic method [46], the 

variational iteration method [47], and the Adomian decomposition method [48]. In addition, 

Khataybeh et al. [42] applied the classical operational matrices of the Bernstein polynomial to 

solve the equation. Also, Parand and Taghavi [49] implemented a collocation method based on a 

rationally scaled generalized Laguerre function to solve the Blasius equation. 

2.3 The Falkner-Skan Equation 

The boundary layer equations are a significant class of nonlinear ordinary differential 

equations with several uses in fluid dynamics and physics [50]. One of these equations is the 

stationary Falkner-Skan boundary layer equation. The Falkner-Skan equation was initially put 

out by Falkner and Skan in 1931[51]. This equation is essential for numerous applications, 

including fluid mechanics, aerospace, heat transfer, glass applications, and polymer 

investigations [20]. 

The third-order ordinary differential equation of the Falkner-Skan equation over a semi-

infinite domain is given by [52]: 

𝑦′′′(𝑥) + 𝑘𝑦(𝑥) 𝑦′′(𝑥) + 𝛽 [𝜖2 − (𝑦′(𝑥))2] = 0,                                                              (6)  

subjected to the boundary conditions: 

𝑦(0) = 0,   𝑦′(0) = 1 − 𝜖,   𝑦′(∞) = 𝜖,                                                                                (7)  

where 𝑘 = 1 is constant. 

The velocity ratio parameter is denoted by 𝜖, and the pressure gradient parameter by 𝛽. The 

Equation (6) is known as the Blasius equation when 𝛽 = 0 and 𝑘 =
1

2
, the Homann flow issue 

when 𝛽 =  
1

2
 and 𝑘 = 1, and the Hiemenz flow problem when 𝛽 =  1 and 𝑘 = 1, see [20]. 

The initial condition 𝑦′′(0) = −0.832666 has been derived by the authors from the boundary 

condition 𝑦′(∞) = 𝜖, using the Padé approximation method [53], and this value will be utilized 

in the current paper. As a result, the following are the initial conditions for the Falkner-Skan 

equation: 

𝑦(0) = 0,   𝑦′(0) = 1 − 𝜖,   𝑦′′(0)  = −0.832666.                                                           (8)  



IHJPAS. 36 (4) 2023 

342 
 

The Falkner-Skan equation has been solved by using a variety of techniques, like the 

homotopy perturbation method [54], the homotopy analysis method [55], the Adomian 

decomposition method [56], the differential transformation method [57], the iterative 

transformation method [58], the Legendre rational polynomials method [59], the shifted 

Chebyshev collocation method [60], and the modified rational Bernoulli functions [61]. 

3. The Basic Concepts of the Proposed Methods 

The fundamental ideas of the suggested methods are presented in this section. In addition, the 

orthogonal polynomials and operational matrices will be introduced as instruments for improving 

the ECM approach to obtain novel approximate solutions to specific nonlinear initial and 

boundary value problems described in section two. 

3.1 The Basic Concepts of the ECM and Their Operational Matrices 

Consider the following 𝑛𝑡ℎ- order ordinary differential equation [29]: 

𝐺(𝑥, 𝑦, 𝑦′, 𝑦′′, … , 𝑦(𝑛)) = 𝑓(𝑥),          𝑎 ≤ 𝑥 ≤ 𝑏,                                                               (9) 

subjected to the initial condition: 𝑦(𝑘)(𝑎) = 𝜇𝑘,       0 ≤ 𝑘 ≤ 𝑛 − 1,                           (10) 

or with the boundary conditions: 𝑦(𝑘)(𝑎) = 𝜔𝑘,   𝑦
(𝑘)(𝑏) = 𝛾𝑘,    0 ≤ 𝑘 ≤

𝑛

2
− 1.        (11)  

Where 𝜇𝑘, 𝜔𝑘, and 𝛾𝑘 are constants and 𝑓(𝑥) is a known function. 

The essential assumption is that the Equation (9) has a unique solution when the initial or 

boundary conditions are determined in the Equations (10) or (11). Moreover, a linear combination 

of 𝑛𝑡ℎ-order functional series based on standard polynomials may be used to represent the 

unknown function 𝑦(𝑥) as follows: 

𝑦(𝑥) = ∑ 𝑎𝑘 𝜓𝑘(𝑥) = 𝜳(𝑥) 𝑨,                                                                                           (12)

𝑛

𝑘=0

 

where 𝜳(𝑥) = [1  𝑥  𝑥2 𝑥3 … 𝑥𝑛] and 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇 , such that 𝑎𝑘, 𝑘 = 0, … , 𝑛, are the 

coefficients, whose values will be specified later. 

Assume 𝜳(𝑥) has the following derivatives: 

𝜳′(𝑥) = 𝜳(𝑥) 𝑩∗ , 𝜳′′(𝑥) = 𝜳(𝑥) (𝑩∗)2, … ,𝜳(𝑛)(𝑥) = 𝜳(𝑥) (𝑩∗)𝑛, 

where 𝑩∗
(𝑛+1)×(𝑛+1), is the operational matrix and its entry values are from the following in the 

standard polynomials: 

𝑩∗ =

[
 
 
 
 
 
0 1 0
0 0 2
0 0 0

⋯
0
0
0

⋮ ⋱ ⋮
0 0 0
0 0 0

⋯
𝑛
0]
 
 
 
 
 

(𝑛+1)×(𝑛+1)

 

Thus, the forms presented below can be used to define the derivatives of the function 𝑦(𝑥): 

𝑦(𝑛)(𝑥) = 𝜳(𝑥) (𝑩∗)𝑛 𝑨,                     where,   𝑛 ≥ 1.                                                       (13) 

Then, the Equations (12) and (13) are substituted into the Equations (9), (10), and (11), to 

provide the following result: 

𝐺(𝑥,   𝜳(𝑥) 𝑨,   𝜳(𝑥) 𝑩∗ 𝑨,   𝜳(𝑥) (𝑩∗)2 𝑨, … ,𝜳(𝑥) (𝑩∗)𝑛 𝑨) = 𝑓(𝑥),                     (14) 

with,  𝜳(𝑎) (𝑩∗)𝑘 𝑨 = 𝜇𝑘 ,     0 ≤ 𝑘 ≤ 𝑛 − 1,                                                                   (15) 

and, 𝜳(𝑎) (𝑩∗)𝑘 𝑨 = 𝜔𝑘,       𝜳(𝑏) (𝑩∗)𝑘 𝑨 = 𝛾𝑘,    0 ≤ 𝑘 ≤
𝑛

2
− 1.                           (16) 

Moreover, the inner product in the Hilbert space 𝐻 = 𝐿2[0,1], is defined as follows: 



IHJPAS. 36 (4) 2023 

343 
 

〈𝑔1, 𝑔2〉 = ∫𝑔1(𝑥) 𝑔2(𝑥) 𝑑𝑥

1

0

.                                                                                                 (17)  

Also, the set of functions 𝜴 = {𝛺0, 𝛺1 … ,𝛺𝑘}, are linearly independent in 𝐻, where 𝛺𝑘 =

𝑥𝑘 , 0 ≤ 𝑘 ≤ 𝑛, is the base function of the standard polynomials [29]. 

Therefore, applying the inner product of the set of base functions 𝜴 with the left and right 

sides of the Equation (14), as given in the Equation (17), yields the matrix equation shown below 

[31]: 

𝑭 = 𝑮,                                                                                                                                          (18) 

where the 𝑖𝑡ℎ row of 𝑭 and 𝑮 in the matrix equation given in the Equation (18) includes the 

following: 

〈𝛺𝑖, 𝐺(𝑥,   𝜳(𝑥) 𝑨,   𝜳(𝑥) 𝑩∗ 𝑨,   𝜳(𝑥) (𝑩∗)2 𝑨,… ,𝜳(𝑥) (𝑩∗)𝑛 𝑨)  〉,   〈𝛺𝑖, 𝑓(𝑥)〉,   0 ≤ 𝑖 ≤ 𝑛. (19)   

Finally, some of the entries in the matrix equation (Equation(18)) will be modified when the 

initial or boundary conditions from the Equations (15) and (16) are substituted. As a result, a 

system of (𝑛 +  1) non-linear algebraic equations are produced, with unknown coefficients 𝑨. 

Then, solve these algebraic equations numerically with applicable programs or sometimes 

analytically. Unique values for the unknown coefficients 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛] can be acquired, 

which are substituted into the Equation (12) to obtain the approximate solution of the Equation 

(9). 

3.2 The Chebyshev Polynomials and Their Operational Matrices 

The following is the definition of the first kind of the Chebyshev polynomials 𝑻𝑛(𝑥) of degree 

𝑛: 

𝑻𝑛(𝑥) = ∑(−1)𝑛−𝑘 2𝑘
(𝑛 + 𝑘 − 1)!

(𝑛 − 𝑘)! (2𝑘)!
(𝑥 + 1)𝑘.                                                         (20)

𝑛

𝑘=0

 

The function 𝑦(𝑥) can be represented by the (𝑛 + 1)-terms of the Chebyshev polynomials of 

the first kind as below [10]: 

𝑦(𝑥) = ∑ 𝑎𝑘 𝑻𝑘(𝑥) = 𝑨𝑇 𝜳(𝑥),                                                                                       (21)

𝑛

𝑘=0

 

where, 𝜳(𝑥) = [𝑻0(𝑥), 𝑻1(𝑥), 𝑻2(𝑥),… , 𝑻𝑛(𝑥)]𝑇 and 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇, such that 𝑎𝑘, 𝑘 =

0, … , 𝑛, are the unknown Chebyshev polynomials coefficients of the first kind, whose values will 

be determined later. 

Moreover, the derivatives of 𝜳(𝑥) can be regarded as: 

𝜳′(𝑥) = 𝑩𝑻
∗  𝜳(𝑥)  ,𝜳′′(𝑥) = (𝑩𝑻

∗ )2 𝜳(𝑥) , … ,𝜳(𝑛)(𝑥) = (𝑩𝑻
∗ )𝑛 𝜳(𝑥), 

where 𝑩𝑻
∗  (𝑛 + 1)× (𝑛 + 1), is the specified derivative's operational matrix, which is defined as 

follows: 

𝑩𝑻
∗ = (𝑑𝑖,𝑗) = {

2𝑖

𝜇𝑗
,       for 𝑗 = 𝑖 − 𝑘,

0             otherwise,

 

where 𝑘 = 1, 3, 5, … , 𝑛 − 1 if 𝑛 is even, or 𝑘 = 1, 3, 5, … , 𝑛 if 𝑛 is odd, 𝜇0 = 2, and 𝜇𝑘 =

1 for all 𝑘 ≥ 1. 

For instance, if 𝑛 is even, the 𝑩𝑻
∗  is written as follows: 



IHJPAS. 36 (4) 2023 

344 
 

𝑩𝑻
∗ =

(

 
 
 
 
 
 

0 0 0 0 0 … 0 0 0
1 0 0 0 0 … 0 0 0
0 4 0 0 0 … 0 0 0
3 0 6 0 0 … 0 0 0
0 8 0 8 0 … 0 0 0
5 0 10 0 10 … 0 0 0
⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 0 0 0

𝑛 − 1 0 2(𝑛 − 1) 0 2(𝑛 − 1) … 2(𝑛 − 1) 0 0
0 2𝑛 0 2𝑛 0 … 0 2𝑛 0)

 
 
 
 
 
 

. 

 

The matrix 𝑩𝑻
∗  is also defined as follows if 𝑛 is odd: 

𝑩𝑻
∗ =

(

 
 
 
 
 

0 0 0 0 … 0 0 0
1 0 0 0 … 0 0 0
0 4 0 0 … 0 0 0
3 0 6 0 … 0 0 0
0 8 0 8 … 0 0 0
⋮ ⋮ ⋮ ⋮ ⋱ 0 0 0
0 2(𝑛 − 1) 0 2(𝑛 − 1) … 2(𝑛 − 1) 0 0
𝑛 0 2𝑛 0 … 0 2𝑛 0)

 
 
 
 
 

. 

Consequently, the derivatives of the function 𝑦(𝑥)  have the following form: 

𝑦(𝑛)(𝑥) = 𝑨𝑇 (𝑩𝑻
∗ )𝑛  𝜳(𝑥),                     where,   𝑛 ≥ 1.                                                   (22) 

3.3 The Bernoulli Polynomials and Their Operational Matrices 

The definition of the Bernoulli polynomials 𝑩𝑛(𝑥) of degree 𝑛 is as follows: 

𝑩𝑛(𝑥) = ∑
𝑛! 𝒃𝑖

  𝑖 !  (𝑛 − 𝑖)! 2𝑛−𝑖
 (𝑥 + 1)𝑛−𝑖

𝑛

𝑖=0

,                                                                     (23) 

where 𝒃𝑖 = 𝑩𝑖(0) is called the Bernoulli number for each 𝑖 = 0,1, …. These numbers are 

calculated by the following identity [62]: 

𝑥

𝑒𝑥 − 1
= ∑𝒃𝑖

∞

𝑖=0

𝑥𝑖

𝑖!
, 

The first few of the Bernoulli numbers are: 𝒃0 = 1, 𝒃1 = −
1

2
, 𝒃2 =

1

6
, 𝒃4 = −

1

30
, …, and 

𝒃2𝑖+1 = 0 for 𝑖 ≥ 1. 

We intend to approximate the solution 𝑦(𝑥) of the problem with the initial or boundary 

conditions in the form [44]: 

𝑦(𝑥) = ∑ 𝑎𝑘 𝑩𝑘(𝑥) =  𝑨𝑇 𝜳(𝑥),                                                                                        (24)

𝑛

𝑘=0

 

where 𝜳(𝑥) = [𝑩0(𝑥),𝑩1(𝑥), 𝑩2(𝑥), … , , 𝑩𝑛(𝑥)]𝑇 and 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇 , such that 𝑎𝑘, 𝑘 =

0, … , 𝑛, are the unknown Bernoulli coefficients, whose values will be identified later. 

Furthermore, the derivatives of 𝜳(𝑥) can be expressed as: 

𝜳′(𝑥) = 𝑩𝓑
∗  𝜳(𝑥) ,𝜳′′(𝑥) = (𝑩𝓑

∗ )2 𝜳(𝑥) , … ,𝜳(𝑛)(𝑥) = (𝑩𝓑
∗ )𝑛 𝜳(𝑥), 

where 𝑩𝓑
∗  (𝑛 + 1)× (𝑛 + 1),  is the differentiation operational Bernoulli matrix, which is defined as 

follows [11]: 



IHJPAS. 36 (4) 2023 

345 
 

𝑩𝓑
∗ =

[
 
 
 
 
0 0 0 . . . 0 0
1 0 0 . . . 0 0
0
⋮
0

2
⋮
0

0
⋮
0

⋯
⋱
. . .

0 0
⋮ ⋮
𝑛 0]

 
 
 
 

(𝑛+1)×(𝑛+1)

 

Then, the derivatives of the function 𝑦(𝑥) can be defined by: 

𝑦(𝑛)(𝑥) = 𝑨𝑇 (𝑩𝓑
∗ )𝑛  𝜳(𝑥),                     where,   𝑛 ≥ 1.                                                   (25) 

3.4 The Laguerre Polynomials and Their Operational Matrices 

The Laguerre polynomials 𝑳𝑛(𝑥) of degree 𝑛 are defined as follows [12]: 

𝑳𝑛(𝑥) = ∑(−1)𝑘
 𝑛!  

(𝑛 − 𝑘)! (𝑘!)2
 𝑥𝑘

𝑛

𝑘=0

 ,                    𝑛 ≥ 0                                                (26) 

with 𝑳𝑛(0) = 1. Moreover, the first few Laguerre polynomials 𝑳𝑛(𝑥) are as follows: 

𝑳0(𝑥) = 1,   𝑳1(𝑥) = 1 − 𝑥,   𝑳2(𝑥) = 1 − 2 𝑥 +
𝑥2

2
,   𝑳3(𝑥) = 1 − 3 𝑥 +

3 𝑥2

2
−

𝑥3

6
,… 

The unknown function 𝑦(𝑥) can be approximated by the (𝑛 + 1)-terms of the Laguerre 

polynomials as [17]: 

𝑦(𝑥) = ∑ 𝑎𝑘 𝑳𝑘(𝑥) = 𝜳(𝑥) 𝑨,                                                                                            (27)

𝑛

𝑘=0

 

where, 𝜳(𝑥) = [𝑳0(𝑥), 𝑳1(𝑥), 𝑳2(𝑥),… , 𝑳𝑛(𝑥)] and 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇 , such that 𝑎𝑘, 𝑘 =

0, … , 𝑛, are the unknown Laguerre coefficients, whose values will be determined later. 

In addition, the relation between the Laguerre polynomials 𝑳𝑛(𝑥) and its integer order 

derivatives is defined by [17]: 

𝜳′(𝑥) = 𝜳(𝑥) 𝑩𝑳
∗  ,𝜳′′(𝑥) = 𝜳(𝑥) (𝑩𝑳

∗)2, … ,𝜳(𝑛)(𝑥) = 𝜳(𝑥) (𝑩𝑳
∗)𝑛, 

where 𝑩𝑳
∗  (𝑛 + 1)× (𝑛 + 1), is the operational matrix of the provided derivative of the Laguerre 

polynomials, which is defined by [17]: 

𝑩𝑳
∗  =

[
 
 
 
 
 
0 −1 −1
0    0 −1
0    0    0

⋯
−1
−1
−1

⋮ ⋱ ⋮
0    0   0
0    0   0

⋯
−1
  0 ]

 
 
 
 
 

(𝑛+1)×(𝑛+1)

. 

Accordingly, the derivatives of the function 𝑦(𝑥) can be expressed by: 

𝑦(𝑛)(𝑥) = 𝜳(𝑥) (𝑩𝑳
∗)𝑛 𝑨,                     where,   𝑛 ≥ 1.                                                       (28) 

4. The Implementation of the ECM and I-ECMs and Numerical Results 

The proposed methods of the ECM and the I-ECMs will be applied in this section to find novel 

approximate solutions, and the numerical results will be presented for three nonlinear problems: 

the Darcy-Brinkman-Forchheimer equation, the Blasius equation, and the Falkner-Skan equation. 

The I-ECMs are based on the base functions of diverse polynomials such as Chebyshev, 

Bernoulli, and Laguerre polynomials, introduced in Equations (20), (23), and (26), respectively, 

with relevant operational matrices. These polynomials are performed in two steps of the proposed 

method's procedures to improve the ECM's accuracy and reliability. First, describe the unknown 

function 𝑦(𝑥) and its derivatives; and second, calculate the inner product to solve the left and 

right sides of the matrix equation explained in Equation (18). 



IHJPAS. 36 (4) 2023 

346 
 

Furthermore, the initial or boundary conditions are substituted, as specified in Equations (15) 

and (16), and some entries of Equation (18) are modified. Therefore, we obtain (𝑛 +  1) 

nonlinear algebraic equations for the unknown coefficients𝑨. By solving this system numerically 

using Mathematica®12, we get the values for the unknown coefficients 𝑎0, 𝑎1, 𝑎2, … , 𝑎𝑛 to obtain 

a novel approximate solution to the nonlinear initial and boundary value problems. 

4.1 Solving the Darcy-Brinkman-Forchheimer Equation by the ECM and I-ECMs 

The ECM and the I-ECMs techniques are used to solve the first problem presented in the 

Equation (1) with boundary conditions in Equation (2). More precisely, for the ECM technique, 

we transform the function 𝑦(𝑥) and its derivatives into matrices by substituting Equations (12) 

and (13) into Equations (1) and (2). Thus, we get the following result: 

𝜳(𝑥) (𝑩∗)2 𝑨 − 𝑠2 (𝜳(𝑥) 𝑨) − 𝐹𝑠(𝜳(𝑥) 𝑨)2 +
1

𝑀
= 0, 

𝜳(0) 𝑩∗ 𝑨 = 0,   𝜳(1) 𝑨 = 0.                                                                                                    (29) 

Then, the procedures have been applied, as shown in Equations (18) and (19), leading to: 

〈𝑥𝑖 ,   𝜳(𝑥) (𝑩∗)2 𝑨 − 𝑠2 (𝜳(𝑥) 𝑨) − 𝐹𝑠(𝜳(𝑥) 𝑨)2〉 = 〈𝑥𝑖 , −
1

𝑀
〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.        (30) 

Substituting Equations (21) and (22) into Equations (1) and (2) for the I-ECMs based on the 

first kind of the Chebyshev polynomials, the following result is obtained: 

𝑨𝑇 (𝑩𝑻
∗ )2  𝜳(𝑥) − 𝑠2 (𝑨𝑇 𝜳(𝑥)) − 𝐹𝑠(𝑨𝑇 𝜳(𝑥))

2
+

1

𝑀
= 0, 

𝑨𝑇 𝑩𝑻
∗   𝜳(0) = 0,   𝑨𝑇 𝜳(1) = 0.                                                                                              (31) 

Additionally, the results of applying Equations (18) and (19) are as follows:  

〈𝑻𝑖(𝑥), 𝑨𝑇 (𝑩𝑻
∗ )2  𝜳(𝑥) − 𝑠2 (𝑨𝑇 𝜳(𝑥)) − 𝐹𝑠(𝑨𝑇 𝜳(𝑥))

2
〉 = 〈𝑻𝑖(𝑥),−

1

𝑀
〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.   (32) 

Implementing the I-ECMs based on the Bernoulli polynomials by substituting Equations (24) 

and (25) into Equations (1) and (2), it follows: 

𝑨𝑇 (𝑩𝓑
∗ )2  𝜳(𝑥) − 𝑠2 (𝑨𝑇 𝜳(𝑥)) − 𝐹𝑠(𝑨𝑇 𝜳(𝑥))

2
+

1

𝑀
= 0, 

𝑨𝑇 𝑩𝓑
∗   𝜳(0) = 0,   𝑨𝑇 𝜳(1) = 0.                                                                                             (33) 

Using the technique described in Equations (18) and (19), the following equation will be given: 

〈𝑩𝑖(𝑥), 𝑨𝑇 (𝑩𝓑
∗ )2  𝜳(𝑥) − 𝑠2 (𝑨𝑇 𝜳(𝑥)) − 𝐹𝑠(𝑨𝑇 𝜳(𝑥))

2
〉 = 〈𝑩𝑖(𝑥),−

1

𝑀
〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.  (34) 

Moreover, applying the I-ECMs based on the Laguerre polynomials by substituting the 

Equations (27) and (28) into the Equations (1) and (2), we obtain: 

𝜳(𝑥) (𝑩𝑳
∗)2 𝑨 − 𝑠2 (𝜳(𝑥) 𝑨) − 𝐹𝑠(𝜳(𝑥) 𝑨)2 +

1

𝑀
= 0, 

𝜳(0) 𝑩𝑳
∗  𝑨 = 0,   𝜳(1) 𝑨 = 0.                                                                                                    (35) 

Subsequently, the procedures as specified in Equations (18) and (19) have been utilized, as 

will be illustrated: 

〈𝑳𝑖(𝑥),𝜳(𝑥) (𝑩𝑳
∗)2 𝑨 − 𝑠2 (𝜳(𝑥) 𝑨) − 𝐹𝑠(𝜳(𝑥) 𝑨)2〉 = 〈𝑳𝑖(𝑥),−

1

𝑀
〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.  (36) 

Additionally, the inner product for the left and right sides of Equations (30), (32), (34), and 

(36), respectively, is used to get the values of 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇 by solving the algebraic 

system of equations. Once the boundary conditions have been applied to Equations (29), (31), 

(33), and (35), respectively, the desired novel approximate solutions are obtained. 



IHJPAS. 36 (4) 2023 

347 
 

If the parameter values are 𝑠 = 1, 𝐹 = 1, and 𝑀 = 1, as in [36], with 𝑛 = 10, then the novel 

approximate solutions for the Darcy-Brinkman-Forchheimer equation will be: 

By applying the ECM based on the standard polynomials: 

𝑦(𝑥) ≈ 0.323852 − 0.285634 𝑥2 + 2.29576 × 10−6𝑥3 − 0.0392379 𝑥4

+ 0.0000786079 𝑥5 + 0.000355229 𝑥6 + 0.000352023 𝑥7

+ 0.0000516179 𝑥8 + 0.00021914 𝑥9 − 0.0000397161 𝑥10. 

Also, by implementing the I-ECMs based on the first kind of the Chebyshev polynomials, we 

obtain: 

𝑦(𝑥) ≈ 0.323852 − 0.285634 𝑥2 + 8.64572 × 10−7𝑥3 − 0.039229 𝑥4 + 0.0000472476 𝑥5

+ 0.000421857 𝑥6 + 0.000264638 𝑥7 + 0.000120859 𝑥8

+ 0.000188736 𝑥9 − 0.0000340335 𝑥10. 

Moreover, by using the I-ECMs based on the Bernoulli polynomials, we achieve: 

𝑦(𝑥) ≈ 0.323852 − 0.285634 𝑥2 + 8.72335 × 10−7𝑥3 − 0.039229 𝑥4 + 0.0000475082 𝑥5

+ 0.000421236 𝑥6 + 0.000265524 𝑥7 + 0.000120109 𝑥8

+ 0.000189083 𝑥9 − 0.0000341013 𝑥10. 

In addition, by utilizing the I-ECMs based on the Laguerre polynomials, we get: 

𝑦(𝑥) ≈ 0.323852 − 0.285634 𝑥2 + 7.16738 × 10−7𝑥3 − 0.0392277 𝑥4

+ 0.0000418568 𝑥5 + 0.000435078 𝑥6 + 0.0002453 𝑥7 + 0.000137559 𝑥8

+ 0.000180872 𝑥9 − 0.0000324758 𝑥10. 

Furthermore, since the exact solution to the Darcy-Brinkman-Forchheimer equation is 

unknown, the 𝑀𝐸𝑅𝑛 has been calculated to determine the accuracy and reliability of the novel 

approximate solution produced by the proposed approaches. The 𝑀𝐸𝑅𝑛 is calculated by: 

𝑀𝐸𝑅𝑛 = 𝑚𝑎𝑥
0≤𝑥≤1

|𝑦′′(𝑥) − 𝑠2𝑦(𝑥) − 𝐹𝑠𝑦2(𝑥) +
1

𝑀
| . 

Figure 2 presents the logarithmic plots for the 𝑀𝐸𝑅𝑛 values obtained by the ECM based on 

the standard polynomials and by the I-ECMs based on the Chebyshev, Bernoulli, and Laguerre 

polynomials, which prove the efficiency and accuracy of these techniques by observation of the 

error values for  𝑛 =  2 to 10, as we found that the error decreases with increasing the values of 

𝑛. 

 

 

 

 

 

Figure 2. Logarithmic plots of 𝑀𝐸𝑅𝑛 to the Darcy-Brinkman-Forchheimer equation. 



IHJPAS. 36 (4) 2023 

348 
 

Also, Figure 3 presents the comparison between the novel approximate solutions calculated 

by the proposed techniques for 𝑛 = 10, 𝑠 = 1, 𝐹 = 1, and 𝑀 = 1. It is evident that impressive 

agreements have been achieved for all the suggested methods. 

 

 

 

 

 

Figure 3. The comparison of the solutions to the Darcy-Brinkman-Forchheimer equation by proposed methods. 

Moreover, the values of the 𝑀𝐸𝑅𝑛 for the novel approximate solutions utilizing ECM and I-

ECMs are also shown in Table 1 with 𝑛 = 10 and parameters 𝑠 = 𝑀 = 1, versus the value of 𝐹, 

which offers the accuracy of these techniques. In addition, it can be observed that the I-ECMs 

based on the Chebyshev polynomial method provide slightly better accuracy with the lowest 

number of errors compared to other techniques. 

 Table 1. The comparison between the 𝑀𝐸𝑅10 when 𝑠 = 𝑀 = 1, and versus the value of 𝐹 for the Darcy-

Brinkman-Forchheimer equation. 

 

4.2 Solving the Blasius Equation by the ECM and I-ECMs 

The ECM and the I-ECMs techniques are utilized to solve the second problem shown in the 

Equations (3) and (5). More precisely, we substitute Equations (12) and (13) into Equations (3) 

and (5) for the technique ECM, converting the function 𝑦(𝑥) and its derivatives as matrices. Thus, 

we obtain the following result: 

𝜳(𝑥) (𝑩∗)3 𝑨 +
1

2
 (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩∗)2 𝑨) = 0, 

𝜳(0) 𝑨 = 0,   𝜳(0) 𝑩∗ 𝑨 = 0,    𝜳(0) (𝑩∗)2 𝑨 = 𝛼.                                                            (37) 

Then, the processes have been applied, as shown in the Equations (18) and (19), so: 

〈𝑥𝑖 ,   𝜳(𝑥) (𝑩∗)3 𝑨 +
1

2
 (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩∗)2 𝑨)〉 = 〈𝑥𝑖 , 0〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.                (38)   

Substituting Equations (21) and (22) into Equations (3) and (5) for the I-ECMs based on the 

first kind of Chebyshev polynomials, it follows: 

𝑨𝑇 (𝑩𝑻
∗ )3  𝜳(𝑥) +

1

2
 (𝑨𝑇  𝜳(𝑥))(𝑨𝑇 (𝑩𝑻

∗ )2  𝜳(𝑥)) = 0, 

𝑨𝑇 𝜳(0) = 0,   𝑨𝑇 𝑩𝑻
∗   𝜳(0) = 0,    𝑨𝑇 (𝑩𝑻

∗ )2  𝜳(0) = 𝛼.                                                  (39) 

And, the results of implementing Equations (18) and (19) are as follows: 

〈𝑻𝑖(𝑥), 𝑨𝑇 (𝑩𝑻
∗ )3  𝜳(𝑥) +

1

2
 (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝑻

∗ )2  𝜳(𝑥))〉 = 〈𝑻𝑖(𝑥), 0〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.          (40)   

𝑭 ECM Standard I-ECMs Chebyshev I-ECMs Bernoulli I-ECMs Laguerre 

2 1.27523 × 10−6 2.75302 × 10−7 2.7852 × 10−7 2.77819 × 10−7 

4 5.33771 × 10−6 1.14476 × 10−6 1.15826 × 10−6 1.17427 × 10−6 

6 0.0000118871 2.52936 × 10−6 2.5595 × 10−6 2.65222 × 10−6 



IHJPAS. 36 (4) 2023 

349 
 

Applying the I-ECMs based on the Bernoulli polynomials by substituting Equations (24) and 

(25) into Equations (3) and (5), the following is obtained: 

𝑨𝑇 (𝑩𝓑
∗ )3  𝜳(𝑥) +

1

2
 (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝓑

∗ )2  𝜳(𝑥)) = 0, 

𝑨𝑇 𝜳(0) = 0,   𝑨𝑇 𝑩𝓑
∗  𝜳(0) = 0,    𝑨𝑇 (𝑩𝓑

∗ )2  𝜳(0) = 𝛼.                                                  (41) 

Using the procedures described in Equations (18) and (19), as a result, the following equation 

will be offered: 

〈𝑩𝑖(𝑥), 𝑨𝑇 (𝑩𝓑
∗ )3  𝜳(𝑥) +

1

2
 (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝓑

∗ )2  𝜳(𝑥))〉 = 〈𝑩𝑖(𝑥), 0〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.        (42) 

Moreover, implementing the I-ECMs based on the Laguerre polynomials by substituting 

Equations (27) and (28) into Equations (3) and (5), it follows that: 

𝜳(𝑥) (𝑩𝑳
∗)3 𝑨 +

1

2
 (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩𝑳

∗)2 𝑨) = 0, 

𝜳(0) 𝑨 = 0,   𝜳(0) 𝑩𝑳
∗  𝑨 = 0,    𝜳(0) (𝑩𝑳

∗)2 𝑨 = 𝛼.                                                            (43) 

Then, the processes have been utilized as given in Equations (18) and (19), which will be 

presented: 

〈𝑳𝑖(𝑥),𝜳(𝑥) (𝑩𝑳
∗)3 𝑨 +

1

2
 (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩𝑳

∗)2 𝑨)〉 = 〈𝑳𝑖(𝑥), 0〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.      (44) 

Furthermore, the inner product for the left and right sides of the Equations (38), (40), (42), and 

(44), respectively, is implemented to obtain the values of 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇by solving the 

algebraic system of equations. Then, the desired novel approximate solutions are achieved by 

applying the initial conditions to Equations (37), (39), (41), and (43), respectively. 

In this problem, we consider the value of 𝛼 = 0.3320573, as in [43] with 𝑛 = 10. The novel 

approximate polynomials for the Blasius equation are: 

By using the ECM based on the standard polynomials: 

𝑦(𝑥) ≈ 0.166029 𝑥2 + 3.40035 × 10−9𝑥3 − 2.07849 × 10−8𝑥4 − 0.000459348 𝑥5

− 1.85524 × 10−7𝑥6 + 2.93847 × 10−7𝑥7 + 2.1901 × 10−6𝑥8 + 2.05083

× 10−7𝑥9 − 8.02077 × 10−8𝑥10. 

Also, by applying the I-ECMs based on the first kind of the Chebyshev polynomials, we 

obtain: 

𝑦(𝑥) ≈ 0.166029 𝑥2 + 2.73849 × 10−10𝑥3 − 4.13564 × 10−9𝑥4 − 0.000459399 𝑥5

− 8.68423 × 10−8𝑥6 + 1.75903 × 10−7𝑥7 + 2.2761 × 10−6𝑥8 + 1.70069

× 10−7𝑥9 − 7.41017 × 10−8𝑥10. 

Moreover, by implementing the I-ECMs based on the Bernoulli polynomials, we achieve: 

𝑦(𝑥) ≈ 0.166029 𝑥2 + 2.44299 × 10−10𝑥3 − 3.8326 × 10−9𝑥4 − 0.000459401 𝑥5

− 8.36908 × 10−8𝑥6 + 1.71538 × 10−7𝑥7 + 2.27964 × 10−6𝑥8

+ 1.68511 × 10−7𝑥9 − 7.38136 × 10−8𝑥10. 

In addition, by utilizing the I-ECMs based on the Laguerre polynomials, we get: 

𝑦(𝑥) ≈ 0.166029 𝑥2 + 1.27153 × 10−10𝑥3 − 2.35939 × 10−9𝑥4 − 0.000459408 𝑥5

− 6.42605 × 10−8𝑥6 + 1.42272 × 10−7𝑥7 + 2.3051 × 10−6𝑥8 + 1.56588

× 10−7𝑥9 − 7.14836 × 10−8𝑥10. 

The exact solution to the Blasius equation is not available. Hence, the 𝑀𝐸𝑅𝑛 has been 

calculated to demonstrate the accuracy of the novel approximate solutions obtained by the 

proposed techniques. The 𝑀𝐸𝑅𝑛 is calculated by: 



IHJPAS. 36 (4) 2023 

350 
 

𝑀𝐸𝑅𝑛 = 𝑚𝑎𝑥
0≤𝑥≤1

|𝑦′′′(𝑥) +
1

2
𝑦(𝑥) 𝑦′′(𝑥)| . 

Figure 4 shows the logarithmic plots for the 𝑀𝐸𝑅𝑛 values obtained by the ECM based on the 

standard polynomials and by the I-ECMs based on the Chebyshev, Bernoulli, and Laguerre 

polynomials, for 𝑛 =  3 to 10, with a value of 𝛼 = 0.3320573, according to previous studies 

[43]. The accuracy and efficiency of these methods can be demonstrated by observing the error 

values for 𝑛, as we observed that the error decreases as the value of 𝑛 increases. 

 

 

 

 

 

 

Figure 4. Logarithmic plots of 𝑀𝐸𝑅𝑛 for the Blasius equation. 

Figure 5 also illustrates the comparison between the novel approximate solutions calculated 

by the proposed techniques for 𝑛 =  10 and 𝛼 = 0.3320573. The figure shows that all of the 

suggested methods have obtained good agreement. 

 

 

 

 

 

Figure 5. The comparison of the solutions for the Blasius equation. 

Table 2 also presents the 𝑀𝐸𝑅𝑛 values for the novel approximate solutions obtained with the 

ECM and the I-ECMs with 𝑛 = 10, illustrating the efficacy of these methods. Furthermore, it can 

be seen that the I-ECMs based on the Laguerre polynomials technique give good accuracy with 

fewer errors compared to the other methods. 

Table 2. The comparison between the 𝑀𝐸𝑅10 for the Blasius equation by proposed methods. 

 

 

𝒏 ECM Standard I-ECMs Chebyshev I-ECMs Bernoulli I-ECMs Laguerre 

10 2.04021 × 10−8 1.64309 × 10−9 1.46579 × 10−9 8.05404 × 10−10 



IHJPAS. 36 (4) 2023 

351 
 

4.3 Solving the Falkner-Skan Equation by the ECM and I-ECMs 

The procedures of the ECM and the I-ECMs techniques can be implemented to solve the third 

problem introduced in Equations (6) and (8). To be more specific, for the ECM technique, we 

transform the unknown function 𝑦(𝑥) with its derivatives as matrices by substituting the 

Equations (12) and (13) into Equations (6) and (8). Thus, we get the following result: 

𝜳(𝑥) (𝑩∗)3 𝑨 + (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩∗)2 𝑨) + 𝛽 [𝜖2 − (𝜳(𝑥) 𝑩∗ 𝑨)2] = 0, 

𝜳(0) 𝑨 = 0,   𝜳(0) 𝑩∗ 𝑨 = 1 − 𝜖,   𝜳(0) (𝑩∗)2 𝑨 = −0.832666.                                  (45) 

Then, the processes have been used as presented in the Equations (18) and (19), so: 

〈𝑥𝑖 , 𝜳(𝑥) (𝑩∗)3 𝑨 + (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩∗)2 𝑨) + 𝛽 [−(𝜳(𝑥) 𝑩∗ 𝑨)2]  〉                 

= 〈𝑥𝑖 , − 𝛽 𝜖2〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.                                                                      (46)   

Substituting Equations (21) and (22) into Equations (6) and (8) for the I-ECMs based on the 

first kind of the Chebyshev polynomials yields the following: 

𝑨𝑇 (𝑩𝑻
∗ )3  𝜳(𝑥) + (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝑻

∗ )2  𝜳(𝑥)) + 𝛽 [𝜖2 − (𝑨𝑇 𝑩𝑻
∗   𝜳(𝑥))2] = 0, 

𝑨𝑇 𝜳(0) = 0,   𝑨𝑇 𝑩𝑻
∗   𝜳(0) = 1 − 𝜖,   𝑨𝑇 (𝑩𝑻

∗ )2  𝜳(0) = −0.832666.                        (47) 

And, by applying the processes shown in the Equations (18) and (19), the following results: 

〈𝑻𝑖(𝑥),  𝑨𝑇 (𝑩𝑻
∗ )3  𝜳(𝑥) + (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝑻

∗ )2  𝜳(𝑥)) + 𝛽 [−(𝑨𝑇 𝑩𝑻
∗   𝜳(𝑥))2]〉

= 〈𝑻𝑖(𝑥),− 𝛽 𝜖2〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.                                                              (48)  

Implementing the I-ECMs based on the Bernoulli polynomials by substituting Equations (24) 

and (25) into Equations (6) and (8), the following is achieved: 

𝑨𝑇 (𝑩𝓑
∗ )3  𝜳(𝑥) + (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝓑

∗ )2  𝜳(𝑥)) + 𝛽 [𝜖2 − (𝑨𝑇 𝑩𝓑
∗   𝜳(𝑥))2] = 0, 

𝑨𝑇 𝜳(0) = 0,   𝑨𝑇 𝑩𝓑
∗   𝜳(0) = 1 − 𝜖,   𝑨𝑇 (𝑩𝓑

∗ )2  𝜳(0) = −0.832666.                        (49) 

Also, by using the techniques as specified in Equations (18) and (19), it follows that: 

〈𝑩𝑖(𝑥),  𝑨𝑇 (𝑩𝓑
∗ )3  𝜳(𝑥) + (𝑨𝑇 𝜳(𝑥))(𝑨𝑇 (𝑩𝓑

∗ )2  𝜳(𝑥)) + 𝛽 [−(𝑨𝑇 𝑩𝓑
∗   𝜳(𝑥))2]〉

= 〈𝑩𝑖(𝑥), − 𝛽 𝜖2〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.                                                                    (50) 

Moreover, applying the I-ECMs based on the Laguerre polynomials by substituting Equations 

(27) and (28) into Equations (6) and (8), we get: 

𝜳(𝑥) (𝑩𝑳
∗)3 𝑨 + (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩𝑳

∗)2 𝑨) + 𝛽 [𝜖2 − (𝜳(𝑥) 𝑩𝑳
∗  𝑨)2] = 0, 

𝜳(0) 𝑨 = 0,   𝜳(0) 𝑩𝑳
∗  𝑨 = 1 − 𝜖,   𝜳(0) (𝑩𝑳

∗)2 𝑨 = −0.832666.                                  (51) 

Then, the processes have been utilized as provided in Equations (18) and (19), which will be 

shown: 

〈𝑳𝑖(𝑥),   𝜳(𝑥) (𝑩𝑳
∗)3 𝑨 + (𝜳(𝑥) 𝑨)(𝜳(𝑥) (𝑩𝑳

∗)2 𝑨) +  𝛽 [−(𝜳(𝑥)  𝑩𝑳
∗  𝑨)2]〉

= 〈𝑳𝑖(𝑥), − 𝛽 𝜖2〉,   ∀ 0 ≤ 𝑖 ≤ 𝑛.                                                                    (52) 

Furthermore, the values of 𝑨 = [𝑎0 𝑎1 𝑎2 …𝑎𝑛]𝑇 are calculated by solving the algebraic 

system of equations obtained by the inner product for the left and right sides of Equations (46), 

(48), (50), and (52), respectively. Then, we utilize the initial conditions to Equations (45), (47), 

(49), and (51), respectively, the desired novel approximate solutions are obtained. 

The novel approximate polynomials for the Falkner-Skan equation when the parameter values 

are as follows: 𝜖 = 0.1, β = 0.5, as in [53], with 𝑛=8, will be: 

By implementing the ECM based on the standard polynomials: 

𝑦(𝑥) ≈ 0.9 𝑥 − 0.416333 𝑥2 + 0.0666511 𝑥3 + 0.0000592155 𝑥4 − 0.00313186 𝑥5

+ 0.000639976 𝑥6 + 0.0000210854 𝑥7 − 0.0000188788 𝑥8. 

Also, by applying the I-ECMs based on the first kind of the Chebyshev polynomials, we 

obtain: 



IHJPAS. 36 (4) 2023 

352 
 

𝑦(𝑥) ≈ 0.9 𝑥 − 0.416333 𝑥2 + 0.0666646 𝑥3 + 0.0000169313 𝑥4 − 0.00305864 𝑥5

+ 0.00056836 𝑥6 + 0.0000581686 𝑥7 − 0.0000267944 𝑥8. 

In addition, by utilizing the I-ECMs based on the Bernoulli polynomials, we get: 

𝑦(𝑥) ≈ 0.9 𝑥 − 0.416333 𝑥2 + 0.0666648 𝑥3 + 0.0000164756 𝑥4 − 0.00305823 𝑥5

+ 0.000568589 𝑥6 + 0.000057627 𝑥7 − 0.0000265718 𝑥8. 

Moreover, by using the I-ECMs based on the Laguerre polynomials, we achieve: 

𝑦(𝑥) ≈ −8.72066 × 10−14 + 0.9 𝑥 − 0.416333 𝑥2 + 0.0666602 𝑥3 + 0.0000532904 𝑥4

− 0.00316698 𝑥5 + 0.00071942 𝑥6 − 0.0000423225 𝑥7 − 9.73011

× 10−7𝑥8. 

Since there is no exact solution to the Falkner-Skan equation, the 𝑀𝐸𝑅𝑛 is computed in order 

to verify the efficiency and accuracy of the novel approximate solutions found by the ECM and 

the I-ECMs. The 𝑀𝐸𝑅𝑛 is calculated by: 

𝑀𝐸𝑅𝑛 = 𝑚𝑎𝑥
0≤𝑥≤1

|𝑦′′′(𝑥) + 𝑦(𝑥) 𝑦′′(𝑥) + 𝛽 [𝜖2 − (𝑦′(𝑥))2]|. 

Figure 6 exhibits the logarithmic plots for the 𝑀𝐸𝑅𝑛 values obtained for the parameters 𝜖 =

0.1, and 𝛽 = 0.5, according to studies [53], by the ECM based on the standard polynomials and 

by the I-ECMs based on the Chebyshev, Bernoulli, and Laguerre polynomials, which demonstrate 

the accuracy and efficiency of these techniques by observing the error values for 𝑛 =  2 to 8. 

We observe that when 𝑛 is increased, the error decreased. 

 

 

 

 

 

 

Figure 6. Logarithmic plots of 𝑀𝐸𝑅𝑛 for the Falkner-Skan equation by proposed methods. 

Also, Table 3 shows the 𝑀𝐸𝑅𝑛 values for the novel approximate solutions achieved with the 

ECM and the I-ECMs with 𝑛 = 8, explaining the accuracy of these methods. Moreover, the I-

ECMs based on the Bernoulli polynomials method, offer slightly better accuracy and fewer errors 

than the other methods. 

 Table 3. The comparison between the 𝑀𝐸𝑅8 for the Falkner-Skan equation by proposed methods. 

 

 

Moreover, Figure 7 shows the comparison between the novel approximate solutions calculated 

by the proposed techniques for 𝑛 =  8, 𝜖 = 0.1, and  𝛽 = 0.5. The figure shows that all of the 

suggested approaches exhibited good agreement. 

𝒏 ECM Standard I-ECMs Chebyshev I-ECMs Bernoulli I-ECMs Laguerre 

8 0.0000935074 0.0000123869 0.0000113839 0.0000385485 



IHJPAS. 36 (4) 2023 

353 
 

  

 

 

 

 

 Figure 7. The comparison of the solutions to the Falkner-Skan equation by proposed methods. 

 

Figure 8. Logarithmic plots of 𝑀𝐸𝑅𝑛 for the Falkner-Skan equation by (a) ECM based on the standard 

polynomials and (b) I-ECMs based on the Chebyshev polynomials. 

 

In addition, Figures (8 and 9) explain the logarithmic plots of the 𝑀𝐸𝑅𝑛 for the novel 

approximate solutions of the Falkner-Skan equation with 𝑛 = 2 to 8, using the ECM and the I-

ECMs when fixed the pressure gradient parameter 𝛽 = 0.5, and increasing the values of the 

velocity ratio parameter as 𝜖 = 0.1, 0.2, 0.3, and 0.4, as chosen in [53]. In Figures (8 and 9), the 

errors decrease when the value of 𝜖 is increased.                                                                                                     

 

 

(a)                                                                                                      (b) 

 

Figure 9. Logarithmic plots of 𝑀𝐸𝑅𝑛 for the Falkner-Skan equation by (a) I-ECMs based on the Bernoulli 

polynomials and (b) I-ECMs based on the Laguerre polynomials. 

 

(a) (b) 



IHJPAS. 36 (4) 2023 

354 
 

Furthermore, Figures (10 and 11) illustrate the logarithmic plots of the 𝑀𝐸𝑅𝑛 for the novel 

approximate solutions of the Falkner-Skan equation with 𝑛 = 2 to 8, by using the ECM and the 

I-ECMs for different values of 𝛽 when fixed the parameter 𝜖 = 0.1. In Figures (10 and 11), it is 

evident that the errors increase as the values of 𝛽 increase. 

 

 

 

 

(a)                                                                                            (b) 

Figure 10. Logarithmic plots of 𝑀𝐸𝑅𝑛 for the Falkner-Skan equation by (a) ECM based on the standard 

polynomials and (b) I-ECMs based on the Chebyshev polynomials. 

 

 

 

 

(a)                                                                                            (b) 

Figure 11. Logarithmic plots of 𝑀𝐸𝑅𝑛 for the Falkner-Skan equation by (a) I-ECMs based on the Bernoulli 

polynomials and (b) I-ECMs based on the Laguerre polynomials. 

 

 

 

5. Conclusions 

In this paper, the effective computational method based on standard polynomials and the novel 

effective computational methods based on the three different types of Chebyshev, Bernoulli, and 

Laguerre polynomials have been implemented to solve three nonlinear models involving initial 

and boundary value problems. Three models, which are well-known nonlinear problems: the 

Darcy-Brinkman-Forchheimer model, the Blasius model, and the Falkner-Skan model, have been 

presented and solved by using our suggested methods. The nonlinear problems are reduced to a 

nonlinear algebraic system of equations solved with Mathematica®12. The novel approximate 



IHJPAS. 36 (4) 2023 

355 
 

solutions were obtained and proved accurate and reliable, even within a few polynomial orders. 

Moreover, the 𝑀𝐸𝑅𝑛 for the proposed methods were calculated. The results show that the 

proposed approaches have higher accuracy and less error. It is also observed that the 𝑀𝐸𝑅𝑛 

results of the proposed methods I-ECMs decrease vastly compared to the ECM. Therefore, the 

proposed novel methods I-ECMs have better accuracy than the ECM. The main conclusion from 

the results is that the Chebyshev polynomials-based I-ECMs have slightly better accuracy than 

the other methods for solving the Darcy-Brinkman-Forchheimer equation. Moreover, the I-ECMs 

based on the Laguerre polynomials are more accurate than the other methods in solving the 

Blasius equation. In addition, the I-ECMs based on the Bernoulli polynomials are slightly more 

accurate than the other methods in solving the Falkner-Skan equation. 

 

References 

1. Hermann, M.; Saravi, M. Nonlinear ordinary differential equations: Analytical 

approximation and numerical methods. Springer India, 2016. 

2. Kounadis, A. N. An efficient and simple approximate technique for solving nonlinear initial 

and boundary-value problems. Computational Mechanics, 1992, 9(3), 221-231. 

3. Talib, I.; Tunc, C.; Noor, Z. A. New operational matrices of orthogonal Legendre 

polynomials and their operational. Journal of Taibah University for Science, 2019, 13(1), 

377-389. 

4. AL-Jawary, M. A.; Salih, O. M. Reliable iterative methods for 1D Swift–Hohenberg 

equation. Arab Journal of Basic and Applied Sciences, 2020, 27(1), 56-66. 

5. Umesh; Kumar, M. Numerical solution of singular boundary value problems using advanced 

Adomian decomposition method. Engineering with Computers, 2021, 37(4), 2853-2863. 

6. Harir, A.; Melliani, S.; El Harfi, H.; Chadli, L. S. Variational iteration method and 

differential transformation method for solving the SEIR epidemic model. International 

Journal of Differential Equations, 2020, 2020, 1-7. 

7. Gohar, M.; Li, C.; Li, Z. Finite difference methods for Caputo–Hadamard fractional 

differential equations. Mediterranean Journal of Mathematics, 2020, 17(6), 1-26. 

8. Roul, P.; Thula, K. A new high-order numerical method for solving singular two-point 

boundary value problems. Journal of Computational and Applied Mathematics, 2018, 343, 

556-574. 

9. Ibraheem, K. I.; Mahmmood, H. S. Algorithm for solving fractional partial differential 

equations using homotopy analysis method with Pade approximation. International Journal 

of Electrical and Computer Engineering, 2022, 12(3), 3335-3342. 

10. Sharma, B.; Kumar, S.; Paswan, M. K.; Mahato, D. Chebyshev operational matrix method 

for Lane-Emden problem. Nonlinear Engineering, 2019, 8(1), 1-9. 

11. Singh, S.; Patel, V. K.; Singh, V. K.; Tohidi, E. Application of Bernoulli matrix method for 

solving two-dimensional hyperbolic telegraph equations with Dirichlet boundary 

conditions. Computers and Mathematics with Applications, 2018, 75(7), 2280-2294. 

12. Bhrawy, A. H.; Taha, T. M. An operational matrix of fractional integration of the Laguerre 

polynomials and its application on a semi-infinite interval. Mathematical Sciences, 

2012, 6(1), 1-7. 

13. AL-Jawary, M. A.; Abdul Nabi, A. J. Three iterative methods for solving Jeffery-Hamel flow 

problem. Kuwait Journal of Science, 2020, 47(1), 1-13.    

14. Agom, E. U.; Ogunfiditimi, F. O.; Bassey, E. V. Homotopy perturbation and Adomian 

decomposition methods on 12th-order boundary value problems. Advances in Mathematics: 

Scientific Journal, 2020, 9(12), 10671-10683. 



IHJPAS. 36 (4) 2023 

356 
 

15. Singh, R. A modified homotopy perturbation method for nonlinear singular Lane–Emden 

equations arising in various physical models. International Journal of Applied and 

Computational Mathematics, 2019, 5(3), 1-15. 

16. 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. 

17. Gürbüz, B.; Sezer, M. Modified operational matrix method for second-order nonlinear 

ordinary differential equations with quadratic and cubic terms. An International Journal of 

Optimization and Control: Theories & Applications, 2020, 10(2), 218-225. 

18. Al-Humedi, H. O.; Al-Saadawi, F. A. The numerical technique based on shifted Jacobi-

Gauss-Lobatto polynomials for solving two dimensional multi-space fractional Bioheat 

equations. Baghdad Science Journal, 2020, 17(4), 1271-1282. 

19. Al-A’asam, J. A. Deriving the composite Simpson rule by using Bernstein polynomials for 

solving Volterra integral equations. Baghdad Science Journal, 2014, 11(3), 1274-1282. 

20. AL-Jawary, M. A.; Ibraheem, G. H. Two meshless methods for solving nonlinear ordinary 

differential equations in engineering and applied sciences. Nonlinear Engineering, 

2020, 9(1), 244-255. 

21. Hasan, P. M.; Sulaiman, N. A. Convergence analysis for the homotopy perturbation method 

for a linear system of mixed Volterra-Fredholm integral equations. Baghdad Science 

Journal, 2020, 17(3 (Suppl.)), 1010-1018. 

22. Ateeah, A. K. Approximate solution for fuzzy differential algebraic equations of fractional 

order using Adomian decomposition method. Ibn AL-Haitham Journal for Pure and Applied 

Science, 2017, 30(2), 202-213. 

23. Abed, S. M.; AL-Jawary, M. A. Efficient iterative methods for solving the SIR epidemic 

model. Iraqi Journal of Science, 2021, 62(2), 613-622. 

24. Sparis, P. D.; Mouroutsos, S. G. The operational matrix of differentiation for orthogonal 

polynomial series. International Journal of Control, 1986, 44(1), 1-15. 

25. Yousefi, S. A.; Behroozifar, M. Operational matrices of Bernstein polynomials and their 

applications. International Journal of Systems Science, 2010, 41(6), 709-716. 

26. Öztürk, Y. Numerical solution of systems of differential equations using operational matrix 

method with Chebyshev polynomials. Journal of Taibah University for Science, 2018, 12(2), 

155-162. 

27. Bazm, S.; Hosseini, A. Bernoulli operational matrix method for the numerical solution of 

nonlinear two-dimensional Volterra–Fredholm integral equations of Hammerstein 

type. Computational and Applied Mathematics, 2020, 39(2), 1-20. 

28. Abdelkawy, M. A.; Taha, T. M. An operational matrix of fractional derivatives of Laguerre 

polynomials. Walailak Journal of Science and Technology, 2014, 11(12), 1041-1055. 

29. Turkyilmazoglu, M. Effective computation of exact and analytic approximate solutions to 

singular nonlinear equations of Lane–Emden–Fowler type. Applied Mathematical 

Modelling, 2013, 37(14-15), 7539-7548.  

30. Turkyilmazoglu, M. An effective approach for numerical solutions of high-order Fredholm 

integro-differential equations. Applied Mathematics and Computation, 2014, 227(15), 384-

398.  

31. Turkyilmazoglu, M. High-order nonlinear Volterra–Fredholm-Hammerstein integro-

differential equations and their effective computation. Applied Mathematics and 

Computation, 2014, 247, 410-416.  

32. Turkyilmazoglu, M. Effective computation of solutions for nonlinear heat transfer problems 

in fins. Journal of Heat Transfer, 2014, 136(9), 091901. 

33. Turkyilmazoglu, M. Solution of initial and boundary value problems by an effective accurate 

method. International Journal of Computational Methods, 2017, 14(06), 1750069.  



IHJPAS. 36 (4) 2023 

357 
 

34. Shaban, M.; Kazem, S.; Rad, J. A. A modification of the homotopy analysis method based 

on Chebyshev operational matrices. Mathematical and Computer Modelling, 2013, 57(5-6), 

1227-1239. 

35. Motsa, S. S.; Sibanda, P.; Shateyi, S. A new spectral-homotopy analysis method for solving 

a nonlinear second order BVP. Communications in Nonlinear Science and Numerical 

Simulation, 2010, 15(9), 2293-2302. 

36. Manafian, J. An optimal Galerkin-homotopy asymptotic method applied to the nonlinear 

second-order BVPs. Proceedings of the Institute of Mathematics and Mechanics, 

2021, 47(1), 156-182. 

37. Waqas, M.; Gulzar, M. M.; Dogonchi, A. S.; Javed, M. A.; Khan, W. A. Darcy–Forchheimer 

stratified flow of viscoelastic nanofluid subjected to convective conditions. Applied 

Nanoscience, 2019, 9(8), 2031-2037. 

38. Awartani, M. M.; Hamdan, M. H. Fully developed flow through a porous channel bounded 

by flat plates. Applied mathematics and computation, 2005, 169(2), 749-757.  

39. Adewumi, A. O.; Akindeinde, S. O.; Aderogba, A. A.; Ogundare, B. S. A hybrid collocation 

method for solving highly nonlinear boundary value problems. Heliyon, 2020, 6(3), 1-10. 

40. Abbasbandy, S.; Shivanian, E.; Hashim, I. Exact analytical solution of forced convection in 

a porous-saturated duct. Communications in Nonlinear Science and Numerical Simulation, 

2011, 16(10), 3981-3989. 

41. El-Gamel, M.; El-Shenawy, A. A numerical solution of Blasius equation on a semi-infinity 

flat plate. SeMA Journal, 2018, 75(3), 475-484. 

42. Khataybeh, S. N.; Hashim, I.; Alshbool, M. Solving directly third-order ODEs using 

operational matrices of Bernstein polynomials method with applications to fluid flow 

equations. Journal of King Saud University-Science, 2019, 31(4), 822-826. 

43. Liao, S. An optimal homotopy-analysis approach for strongly nonlinear differential 

equations. Communications in Nonlinear Science and Numerical Simulation, 2010, 15(8), 

2003-2016. 

44. Rani, D.; Mishra, V. Numerical inverse Laplace transform based on Bernoulli polynomials 

operational matrix for solving nonlinear differential equations. Results in Physics, 2020, 16, 

102836. 

45. Yao, B.; Chen, J. A new analytical solution branch for the Blasius equation with a shrinking 

sheet. Applied Mathematics and Computation, 2009, 215(3), 1146-1153. 

46. Marinca, V.; Herişanu, N. The optimal homotopy asymptotic method for solving Blasius 

equation. Applied Mathematics and Computation, 2014, 231, 134-139. 

47. Wazwaz, A. M. The variational iteration method for solving two forms of Blasius equation 

on a half-infinite domain. Applied Mathematics and Computation, 2007, 188(1), 485-491. 

48. Abbasbandy, S. A numerical solution of Blasius equation by Adomian’s decomposition 

method and comparison with homotopy perturbation method. Chaos, Solitons & Fractals, 

2007, 31(1), 257-260. 

49. Parand, K.; Taghavi, A. Rational scaled generalized Laguerre function collocation method 

for solving the Blasius equation. Journal of Computational and Applied Mathematics, 

2009, 233(4), 980-989. 

50. Abbasbandy, S.; Hajishafieiha, J. Numerical solution to the Falkner-Skan equation: a novel 

numerical approach through the new rational a-polynomials. Applied Mathematics and 

Mechanics, 2021, 42(10), 1449-1460. 

51. Falkner, V. M.; Skan, S. W. Solutions of the boundary-layer equations. The London, 

Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 1931, 12(80), 865-

896. 

52. Bakodah, H. O.; Ebaid, A.; Wazwaz, A. M. Analytical and numerical treatment of Falkner-

Skan equation via a transformation and Adomian’s method. Romanian Reports in Physics, 

2018, 70(2), 1-17. 



IHJPAS. 36 (4) 2023 

358 
 

53. AL-Jawary, M. A.; Adwan, M. I. Reliable iterative methods for solving the Falkner-Skan 

equation. Gazi University Journal of Science, 2020, 33(1), 168-186. 

54. Moallemi, N.; Shafieenejad, I.; Hashemi, S. F.; Fata, A. Approximate explicit solution of 

Falkner-Skan equation by homotopy perturbation method. Research Journal of Applied 

Sciences, Engineering and Technology, 2012, 4(17), 2893-2897. 

55. Abbasbandy, S.; Hayat, T. Solution of the MHD Falkner-Skan flow by homotopy analysis 

method. Communications in Nonlinear Science and Numerical Simulation, 2009, 14(9-10), 

3591-3598. 

56. Elgazery, N. S. Numerical solution for the Falkner-Skan equation. Chaos, Solitons & 

Fractals, 2008, 35(4), 738-746. 

57. Kuo, B. L. Application of the differential transformation method to the solutions of Falkner-

Skan wedge flow. Acta Mechanica, 2003, 164(3), 161-174. 

58. Temimi, H.; Ben-Romdhane, M. Numerical solution of Falkner-Skan equation by iterative 

transformation method. Mathematical Modelling and Analysis, 2018, 23(1), 139-151. 

59. Guo, B. Y.; Shen, J.; Wang, Z. Q. A rational approximation and its applications to differential 

equations on the half line. Journal of scientific computing, 2000, 15(2), 117-147. 

60. Kajani, M. T.; Maleki, M.; Allame, M. A numerical solution of Falkner-Skan equation via a 

shifted Chebyshev collocation method. AIP Conference Proceedings, 2014, 1629(1), 381-

386. 

61. Calvert, V.; Razzaghi, M. Solutions of the Blasius and MHD Falkner-Skan boundary-layer 

equations by modified rational Bernoulli functions. International Journal of Numerical 

Methods for Heat & Fluid Flow, 2017, 27(8), 1687-1705. 

62. Azodi, H. D.; Yaghouti, M. R. Bernoulli polynomials collocation for weakly singular 

Volterra integro-differential equations of fractional order. Filomat, 2018, 32(10), 3623-

3635.