358 

 IHJPAS. 37 (1) 2024 

Ibn Al-Haitham Journal for Pure and Applied Sciences 

Journal homepage: jih.uobaghdad.edu.iq 

PISSN: 1609-4042, EISSN: 2521-3407 

 

 

 

                 

1Amna M. Mahdi*     2Majeed A. AL-Jawary     3Mustafa Turkyilmazoglu     

University of Haitham, -r Pure Sciences Ibn ALDepartment of Mathematics, College of Education fo 1,2

Baghdad, Baghdad, Iraq. 

3Department of Mathematics, Hacettepe University, Ankara, Turkey. 

China Medical University, Department of Medical Research, China Medical University Hospital, 3

.Taichung, Taiwan 

 
*Corresponding Author: amna.meqdad1203a@ihcoedu.uobaghdad.edu.iq 

 

 

 

 

Abstract 

The method of operational matrices based on different types of polynomials such as Bernstein, 

shifted Legendre and Bernoulli polynomials is introduced and implemented to solve the nonlinear 

Blasius equations approximately. The nonlinear differential equation is converted into a system of 

nonlinear algebraic equations that can be solved using Mathematica®12. The efficiency of these 

methods has been studied by calculating the maximum error remainder (,𝑀𝐸𝑅-𝑛.), and it was 

found that their efficiency increases with increasing polynomial degree (n) as the errors decrease. 

Moreover, the approximate solutions obtained by the proposed methods are compared with the 

solution of the fourth-order Runge-Kutta method (RK4), which gives a very good agreement. In 

addition, the convergence of the proposed approximation methods is given based on one of the 

results of the Banach fixed point theorem. 

 

Keywords: Blasius equations; Bernstein polynomial; Legendre polynomial; Bernoulli 

polynomial; operational matrices. 

 

1. Introduction 

         The nonlinear ordinary differential equations (NODE) have many practical applications in 

engineering and applied sciences, such as fluid flow, current flow in electric circuits, heat 

dissipation in solid bodies, seismic wave propagation, population increase or decrease, and many 

other topics [1,2]. In approximation theory and numerical analysis, polynomials are particularly 

useful tools [3]. Consequently, the polynomial series then the operational matrices (OM) are used 

to simplify the unknown function by transforming it into a system of algebraic equations that can 

be easily solved without integration and differentiation. Several studies have been conducted on 

this technique, which uses OM methods to solve many problems based on various polynomials.  

Many researchers have solved various problems using OM based on the Bernstein polynomial 

(BOM), such as: [4] solved odd boundary value problems, [5] studied third order equations ODE, 

[6] solved fractional integral equations. Moreover, there are many researchers who have used 

operational matrices based on different polynomials, such as [7] used the Legendre operational 

Novel Approximate Solutions for Nonlinear Blasius Equations 

Received 24 February 2023, Received 2  May 2023, Accepted 8 May 2023, Published 20 January 2024 

 

 
doi.org/10.30526/37.1.3292 

https://jih.uobaghdad.edu.iq/index.php/j/index#1609-4042
https://jih.uobaghdad.edu.iq/index.php/j/index#2521-3407
mailto:amna.meqdad1203a@ihcoedu.uobaghdad.edu.iq
https://orcid.org/0000-0001-7140-2144
mailto:amna.meqdad1203a@ihcoedu.uobaghdad.edu.iq
https://orcid.org/0000-0003-3967-0012
mailto:majeed.a.w@ihcoedu.uobaghdad.edu.iq
https://orcid.org/0000-0003-0412-4580
mailto:turkyilm@hacettepe.edu.tr


IHJPAS. 37 (1) 2024 

359 
 

matrix (LOM) method to solve the fractional-order two-dimensional integral equations. In [8] 

Sharma et al. solved the Lane-Emden equations using the Chebyshev operational matrix (ChOM) 

method. Also, OM with Bernoulli polynomials (BrOM) was used by Bazm [9] to solve some types 

of integral equations. [10] studied the magnetohydrodynamic squeezing fluid, straight fin problem, 

Jeffery-Hamel flow, and Falkner-Skan equation using BOM and ChOM. They also studied the 

solution of nonlinear thin- film flow of 3rd-grade fluid problems with LOM [11]. In [12-20], other 

types of polynomials were used to solve different types of problems. 

          There is an application for third-order NODEs that occurs in fluid mechanics with the 

laminar viscous flow and the various aspects of the hydrodynamic boundary layer problem is 

the nonlinear Blasius equation [21,22]. Several researchers have worked on the solution of this 

equation: [23] obtained a numerical solution by the Runge-Kutta method. [21] solved the 

Blasius equation, the Duffing equation, the Van der Pol equation, and the Jerk equation using 

the numerical approach of the inverse Laplace transform based on the BrOM integration 

technique. Also, He [24] used the variational iteration method to obtain an analytic approximate 

solution for the Blasius equation. [25] solved it using Adomian’s decomposition method. 

Moreover, [26] studied it with the new homotopy perturbation method.  

The aim of this paper is to use the method of operational matrices based on different types of 

polynomials such as Bernstein, shifted Legendre and Bernoulli polynomials will be used to 

solve the nonlinear Blasius equations and novel approximate solutions will be obtained.  

The structure of this paper is as follows: In Section two, the Blasius equation is introduced. 

Section three gives the orthogonal polynomials on which the operational matrices depend, 

namely the Bernstein polynomial, the shifted Legendre polynomial, and the Bernoulli 

polynomial. In Section four, the convergence of the proposed methods is explained. In section 

five, the problem will be solved using the proposed methods. Finally, Section six presents the 

conclusions. 

2. The Blasius Equation 

This equation is important because it appears in many hydrodynamic boundary layer problems 

as well as in the fluid mechanics of laminar viscous flows and its formula is given by [22]: 

 𝑓′′′(𝑥) +
1

2
𝑓(𝑥)𝑓′′(𝑥) = 0,                                                                                                     (1) 

with boundary conditions: 𝑓(0) = 𝑓′(0) = 0, 𝑓′(∞) = 1.                                                       (2) 

To solve this equation, the boundary conditions were converted to initial conditions by 

computing  𝑓′′(0) instead of  𝑓′(∞), as Liao [27] did and found 𝑓′′(0) = 0.3320573, followed 

by Khataybeh et al. [5] used the value 𝑓′′(0) = 1. 

In general, the Blasius equation represents a model of the attitude of a two-dimensional stable 

laminar viscous flow on a semi-infinite flat plate in which the flowing fluid is incompressible, 

but the boundary layer assumption must be governed by the continuity and Navier-Stokes 

equations of motion, but it should be noted that the fluid flow velocity decreases sharply from 

𝑈 to 0 at 𝑦 = 0 as 𝑥 changes from −0 to + 0. 

 

 



IHJPAS. 37 (1) 2024 

360 
 

3. The operational matrices of the orthogonal polynomials 

Orthogonal polynomials play an important role in pure and applied mathematics as well as in 

numerical computation [28]. Three types of these polynomials will be used: Bernstein, shifted 

Legendre, and Bernoulli to solve the Blasius equation. 

3.1   The Bernstein polynomials 

The Bernstein Polynomials of 𝑛𝑡ℎ degree in [0,1] are defined by [4,10,29]: 

𝐵𝑖,𝑛(𝑥) = (
𝑛
𝑖
) 𝑥𝑖  (1 − 𝑥)𝑛−𝑖,   𝑖 = 0,1,2, … , 𝑛,                                                                      (3) 

where (
𝑛
𝑖
) =

𝑛!

𝑖! (𝑛−𝑖)!
. 

Or the recursive definition over [0,1] is used to generate these polynomials, resulting in 

Bernstein polynomial being represented as follows: 

𝐵𝑖,𝑛(𝑥) = (1 − 𝑥) 𝐵𝑖,𝑛−1(𝑥) + 𝑥 𝐵𝑖−1,𝑛−1(𝑥).                                                                        (4) 

Practically only the first (𝑛 + 1) terms of the Bernstein polynomials of degree 𝑛 are satisfied, 

because 𝐵𝑖,𝑛(𝑥) = 0 if  𝑖 < 0 or 𝑛 < 𝑖. There are many properties make Bernstein polynomial 

important, some of them are: 

i. Property of Positivity: 𝐵𝑖,𝑛(𝑥) > 0 for all 0 < 𝑖 < 𝑛 and all 𝑥𝜖[0,1]. 

ii. Unity partition property: ∑ 𝐵𝑖,𝑛(𝑥) = ∑ 𝐵𝑖,𝑛−1(𝑥) =. . . = ∑ 𝐵𝑖,1(𝑥) = 1.
1
𝑖=0

𝑛−1
𝑖=0

𝑛
𝑖=0  

Moreover, the type of linear combination shown below can be used to approximate any 

polynomial of 𝑛𝑡ℎdegree in 𝐿2[0,1]: 

𝑓(𝑥) = ∑ 𝑐𝑖 𝐵𝑖,𝑛(𝑥)
𝑛
𝑖=0 = 𝐶𝑇𝜙(𝑥),                                                                                          (5) 

 where 𝐶𝑇 = [𝑐0, 𝑐1, 𝑐2, … , 𝑐𝑛], and 𝜙(𝑥) = [𝐵0,𝑛, 𝐵1,𝑛, 𝐵2,𝑛, … , , 𝐵𝑛,𝑛]
𝑇
. 

In addition, 𝜙(𝑥) can be decomposed as the product of a (𝑛 + 1) × (𝑛 + 1) matrix 𝐴 and a 

(𝑛 + 1) × 1 vector 𝑋, i.e.:  

𝜙(𝑥) = 𝐴 𝑋,                                                                                                                              (6) 

where: 

𝐴 =

[
 
 
 
 

 

(−1)0 (𝑛
0
) (−1)1 (𝑛

0
) (𝑛−0

1
) . . . (−1)𝑛−0 (𝑛

0
) (𝑛−0

𝑛−0
)

0 (−1)0 (𝑛
𝑖
) . . . (−1)𝑛−𝑖 (𝑛

𝑖
) (𝑛−𝑖

𝑛−𝑖
)

⋮
0

⋮
0

⋱
. . .
                  

⋮
(−1)0 (𝑛

𝑛
) ]
 
 
 
 

 , X=

[
 
 
 
 

 

1
𝑥
𝑥2 
⋮
𝑥𝑛 ]
 
 
 
 

.                                         

(7) 

The determinant of the matrix 𝐴 is |𝐴| = ∏ (𝑛
𝑖
)𝑛

𝑖=0 . As a result, 𝐴 is an invertible matrix. 

Let us introduce the BOM. If 𝐷𝐵 is the OM of derivative of size (𝑛 + 1) × (𝑛 + 1), then the 

derivative of 𝜙(𝑥) is: 

𝑑𝜙(𝑥)

𝑑𝑥
= 𝐷𝐵𝜙(𝑥);  𝑥𝜖[0,1],                                                                                                         (8) 



IHJPAS. 37 (1) 2024 

361 
 

where 𝐷𝐵 is given as 𝐷𝐵 = 𝐴 𝑉 𝐵
∗ ,where: 

𝑉 =

[
 
 
 
 
0 0 0 ⋯ 0
1 0 0 ⋯ 0
0
⋮
0

2
⋮
0

0
⋮
0

⋯
⋱
⋯

0
⋮
𝑛]
 
 
 
 

(𝑛+1)×𝑛

,𝐵∗ =

[
 
 
 
 

 

𝐴1
−1

𝐴2
−1

𝐴3
−1

⋮
𝐴𝑛
−1

 

]
 
 
 
 

𝑛×(𝑛+1)

,                                                           (9) 

where 𝐵∗ represents the 𝛼𝑡ℎ rows of 𝐴−1for 𝛼=1, 2, …, n. 

Furthermore, the generalization of (8) can be written as follows: 

𝑑𝑛𝜙(𝑥)

𝑑𝑥𝑛
= (𝐷𝐵)

𝑛 𝜙(𝑥);  𝑛 = 1, 2, . ..  .                                                                                      (10) 

Thus, it will be able to approximate the derivatives of 𝑓(𝑥) in terms of OM as follows: 

𝑓′(𝑥) = 𝐶𝑇𝐷𝐵  𝜙(𝑥), 𝑓′′(𝑥) = 𝐶
𝑇(𝐷𝐵)

2 𝜙(𝑥), . . ., 𝑓(𝑛)(𝑥) = 𝐶𝑇(𝐷𝐵)
𝑛 𝜙(𝑥).                  (11) 

This approximation and Eq.(5) are also applied to all conditions of the equation. 

To solve the equation, we replace 𝑓(𝑥) and its derivatives by Eqs.(5) and (11), then modify 𝑥 

by appropriate points from Chebyshev roots, called collocation nodes, as follows: 

𝑥𝑟 =
1

2
 (𝑐𝑜𝑠

𝑟 𝜋

𝑛
+ 1),   𝑟 = 1, . . . , 𝑛 − 1.                                                                                (12) 

This results in a system of algebraic equations that can be solved using software such as 

Mathematica or MATLAB to obtain the value of the vector 𝐶𝑇 in Eq.(5). 

3.2   The shifted Legendre polynomials 

   The Legendre polynomials of 𝑛𝑡ℎ order are defined on [-1,1] as [11,13]: 

𝐿0(𝑡) = 1, 𝐿1(𝑡) = 𝑡, . . ., 𝐿𝑛+1(𝑡) =
2 𝑛+1

𝑛+1
 𝑡 𝐿𝑛(𝑡) −

𝑛

𝑛+1
 𝐿𝑛−1(𝑡), 𝑛 = 1,2, . ..   .           (13)  

Legendre polynomials can be defined on [0,1] by changing the variable 𝑡 = 2𝑥 − 1 and 

denoting it as a shifted Legendre polynomial denoted by 𝑃𝑛(𝑥), and then calculated as follows: 

𝑃0(𝑥) = 1,  𝑃1(𝑥) = 2𝑥 − 1, . . ., 𝑃𝑛+1(𝑥) =
(2𝑛+1)(2𝑥−1)

𝑛+1
𝑃𝑛(𝑥) −

𝑛

𝑛+1
𝑃𝑛−1(𝑥), 𝑛 = 1,2, . ..  .  

(14)  

Also, 𝑃𝑛(𝑥) can be written as shown:  

𝑃𝑛(𝑥) = ∑
(−1)𝑛+𝑖 (𝑛+𝑖)

(𝑛−𝑖)! (𝑖!)2
 𝑥𝑖 .𝑛

𝑖=0                                                                                                  (15) 

Any function 𝑓(𝑥) ∈ 𝐿2[0,1], could be characterized by a shifted Legendre polynomial as 

follows: 

𝑓(𝑥) = ∑ 𝑐𝑖 𝑃𝑖(𝑥)
∞
𝑖=0 ,                                                                                                              (16) 

where 𝑐𝑖 = (2 𝑖 + 1) ∫ 𝑓(𝑥)
1

0
𝑃𝑖(𝑥) 𝑑𝑥;  𝑖 = 0, 1, 2, . ..  .  

In practice, only the first (𝑛 + 1) terms of Eq.(16) will be considered: 

𝑓(𝑥) = ∑ 𝑐𝑖 𝑃𝑖(𝑥)
𝑛
𝑖=0 = 𝐶𝑇𝜙(𝑥),                                                                                              (17) 



IHJPAS. 37 (1) 2024 

362 
 

where 𝐶𝑇 = [𝑐0, 𝑐1, . . . , 𝑐𝑛], and 𝜙(𝑥) = [𝑃0(𝑥), 𝑃1(𝑥), . . . , 𝑃𝑛(𝑥)]
𝑇 . 

The derivatives of 𝜙(𝑥)can be defined as: 

𝑑𝜙(𝑥)

𝑑𝑥
= 𝐷𝐿𝜙(𝑥),

𝑑2𝜙(𝑥)

𝑑𝑥2
= (𝐷𝐿)

2𝜙(𝑥), . . .,
𝑑𝑛𝜙(𝑥)

𝑑𝑥𝑛
= (𝐷𝐿)

𝑛𝜙(𝑥                                        (18) 

where 𝐷𝐿 is the (𝑛 +  1) × (𝑛 +  1) OM of derivative, given as: 

𝐷𝐿 = (𝑑𝑖𝑗) {
2 (2𝑗 + 1),   𝑗 = 𝑖 − 𝑘,   {

𝑘 = 1,3, . . . , 𝑛,   if 𝑛 odd,
𝑘 = 1,3, . . . , 𝑛 − 1,   if 𝑛 even,

0                                                otherwise.
                                  (19) 

Then write the derivatives of 𝑓(𝑥) with respect to the OM as follows: 

𝑓′(𝑥) = 𝐶𝑇𝐷𝐿𝜙(𝑥), 𝑓′′(𝑥) = 𝐶
𝑇(𝐷𝐿)

2 𝜙(𝑥), . . ., 𝑓(𝑛)(𝑥) = 𝐶𝑇(𝐷𝐿)
𝑛 𝜙(𝑥).                    (20) 

This helps to solve the NODE by substituting 𝑓(𝑥) and its derivatives as in Eqs.(17) and (20). 

Also, compensated in the conditions with the NODE. Then the collocation nodes are inserted 

into these equations to produce a system of algebraic equations (𝑛 +  1) that can be solved 

with software such as MATLAB or Mathematica to obtain the coefficients of the vector 𝐶𝑇 . 

3.3   The Bernoulli polynomials 

The 𝑛𝑡ℎBernoulli polynomials on [0,1] are defined as follows [9,21]: 

𝐵𝑟𝑛(𝑥) = ∑ (𝑛
𝑖
)𝐵𝑟𝑖 𝑥

𝑛−𝑖𝑛
𝑖=0 ,                                                                                                   (21) 

where 𝐵𝑟𝑖 = 𝐵𝑟𝑖(0) is the Bernoulli number for all 𝑖 = 0,1, . . . , 𝑛. These numbers are 

calculated by: 

       𝐵𝑟𝑛 = −∑
(−1)𝑖

𝑖

𝑛+1
𝑖=1 (𝑛+1

𝑖
)∑ 𝑗𝑛𝑖

𝑗=1 ,                                                                                   (22) 

For 𝑛 ≥ 0, 𝑛 ≠ 1. If 𝑛 = 1, then 𝐵𝑟1 = −
1

2
. 

Therefore, some Bernoulli polynomials can be represented as follows: 

𝐵𝑟0(𝑥) = 1,    𝐵𝑟1(𝑥) = 𝑥 −
1

2
,    𝐵𝑟2(𝑥) = 𝑥

2 − 𝑥 +
1

6
,    𝐵𝑟3(𝑥) = 𝑥

3 −
3

2
𝑥2 +

1

2
x,  

𝐵𝑟4(𝑥) = 𝑥
4 − 2𝑥3 + 𝑥2 −

1

30
.  

Bernoulli polynomials are used in many fields of mathematics, which led to the discovery of 

important properties for them, some of which we mention below: 

i. 
𝑑𝐵𝑟𝑛(𝑥)

𝑑𝑥
= 𝑛 𝐵𝑟𝑛−1(𝑥), 𝑛 ≥ 1. 

ii. ∫ 𝐵𝑟𝑛(𝑥) 𝑑𝑥 =
𝐵𝑟𝑛+1(𝑧)−𝐵𝑟𝑛+1(𝑎)

𝑛+1

𝑧

𝑎
. 

iii. ∫ 𝐵𝑟𝑛(𝑥) 𝑑𝑥 = 0,
1

0
𝑛 ≥ 1. 

iv. ∫ 𝐵𝑟𝑛(𝑥) 𝐵𝑟𝑚(𝑥) 𝑑𝑥 = (−1)
𝑛−11

0

𝑛! 𝑚!

(𝑛+𝑚)!
 𝐵𝑟𝑛+𝑚. 

v. 𝐵𝑟𝑛(1 − 𝑥) = (−1)
𝑛𝐵𝑟𝑛(𝑥). 

vi. 𝐵𝑟(𝑥 + 1)−𝐵𝑟𝑛(𝑥) = 𝑛 𝑥
𝑛−1. 



IHJPAS. 37 (1) 2024 

363 
 

Therefore, it is now easy to express any 𝑓(𝑥) ∈ 𝐿2[0,1] by the linear combination of the 

Bernoulli polynomial as: 

𝑓(𝑥) = ∑ 𝑐𝑖
𝑛
𝑖=0 𝐵𝑟𝑖(𝑥) = 𝐶

𝑇  𝐵𝑟(𝑥),                                                                                       (23) 

where 𝐶𝑇 = [𝑐0, 𝑐1, … , 𝑐𝑛], and 𝐵𝑟(𝑥) = [ 𝐵𝑟0(𝑥), 𝐵𝑟1(𝑥), . . . , 𝐵𝑟𝑛(𝑥)]
𝑇. For all 𝑖 = 0,1, … , 𝑛, 

the matrix form of 𝐵𝑟𝑖(𝑥) can be obtained as follows: 

𝐵𝑟𝑖(𝑥) =∑(
𝑖

𝑗
)𝐵𝑟𝑗  𝑥

𝑖−𝑗

𝑖

𝑗=0

, 

             = (
𝑖

𝑖
) 𝐵𝑟𝑖 𝑥

0 + (
𝑖

𝑖 − 1
)𝐵𝑟𝑖−1 𝑥 + ⋯+ (

𝑖

1
)𝐵𝑟1 𝑥

𝑖−1 + (
𝑖

0
)𝐵𝑟0 𝑥

𝑖, 

= [(
𝑖

𝑖
) 𝐵𝑟𝑖    (

𝑖

𝑖 − 1
)𝐵𝑟𝑖−1 ⋯    (

𝑖

1
)𝐵𝑟1 (

𝑖

0
)𝐵𝑟0      0   0   . . .   0⏞      

𝑛−𝑖

]  

[
 
 
 
 
 
1…
⋮
 𝑥𝑖

 𝑥𝑖+1

⋮
 𝑥𝑛 ]

 
 
 
 
 

,       

           = 𝑀𝑖 𝐻(𝑥).                                                                                                                     (24) 

where 

𝐻(𝑥) = [1, 𝑥, . . . , 𝑥𝑛]𝑇and 𝑀𝑖 = [(𝑖
𝑖
)𝐵𝑟𝑖   (

𝑖
𝑖−1
)𝐵𝑟𝑖−1 ⋯    (𝑖

1
)𝐵𝑟1 (𝑖

0
)𝐵𝑟0      0   0   . . .    0⏞      

𝑛−𝑖

]. 

As a result, 𝐵𝑟(𝑥) in Eq.(23) can be written as follows: 

𝐵𝑟(𝑥) = [𝐵𝑟0(𝑥), 𝐵𝑟1(𝑥),… , 𝐵𝑟𝑛(𝑥)]
𝑇 ,  

            = [𝑀0𝐻(𝑥),𝑀1𝐻(𝑥),… ,𝑀𝑛𝐻(𝑥)]
𝑇,   

              = [𝑀0, 𝑀1, … ,𝑀𝑛]
𝑇𝐻(𝑥),  

            =

[
 
 
 
 
(0
0
)𝐵𝑟0             0                   0                 ⋯ 0

(1
1
)𝐵𝑟1        (1

0
)𝐵𝑟0             0                  ⋯ 0

⋮   
(𝑛
𝑛
)𝐵𝑟𝑛

 
⋮

    ( 𝑛
𝑛−1
)𝐵𝑟𝑛−1

⋮             
( 𝑛
𝑛−2
)𝐵𝑟𝑛−2

⋯
⋯

⋮
(𝑛
0
)𝐵𝑟0]

 
 
 
 

 𝐻(𝑥),  

            = 𝑀 𝐻(𝑥).                                                                                                              (25) 

While the derivative of it will be: 

       
𝑑𝐵𝑟(𝑥)

𝑑𝑥
= 𝐷𝐵𝑟 𝐵𝑟(𝑥),                                                                                                    (26) 

where 𝐷𝐵𝑟 =

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

2
⋮
0

0
⋮
0

⋯
⋱
. . .

0 0
⋮ ⋮
𝑛 0]

 
 
 
 

 is the (𝑛 + 1)  ×  (𝑛 + 1) OM. 

 

 



IHJPAS. 37 (1) 2024 

364 
 

Moreover, the 𝑛𝑡ℎderivative: 

     
𝑑𝑛𝐵𝑟(𝑥)

𝑑𝑥𝑛
= (𝐷𝐵𝑟)

𝑛 𝐵𝑟(𝑥).                                                                                                   (27) 

To complete the solution, we must find the value of the vector 𝐶𝑇 in Eq.(23), using the same 

procedures as the previous methods. 

 

4. The convergence of the proposed methods 

In this section, the convergence analysis for the proposed methods and fundamental theorem 

are discussed. 

Theorem 4.1.   Let a Banach space 𝐴 ⊂ 𝑅 be given with a norm ‖ ‖ defined on it. Taking 𝑓1(𝑥) 

as an approximate solution obtain from the first iteration n, we construct the following sequence 

regarding the solution of Eq.(1): 

𝑣1(𝑥) = 𝑓1(𝑥), 𝑣𝑖(𝑥) = 𝑓𝑖(𝑥) − 𝑓𝑖−1(𝑥), (𝑖 ≥ 2). 

Then, the assertions are: 

(i) Provided that for all 𝑖 there exist 0 < 𝛽𝑖 < 1 such that ‖𝑣𝑖+1(𝑥)‖ ≤ 𝛽𝑖‖𝑣𝑖(𝑥)‖, the 

series ∑ 𝑣𝑖(𝑥)
∞
𝑖=1  is than convergent and so 𝑓(𝑥) = ∑ 𝑣𝑖(𝑥)

∞
𝑖=1  in the interval of 

interest 𝑥 ∈ Γ. 

(ii) Otherwise, for all 𝑖 there exist 𝛽𝑖 > 1 leading to ‖𝑣𝑖+1(𝑥)‖ ≥ 𝛽𝑖‖𝑣𝑖(𝑥)‖, the series 

∑ 𝑣𝑖(𝑥)
∞
𝑖=1  and thus, the proposed method diverges in the interval of interest 𝑥 ∈ Γ. 

Proof: See [30]. 

 

Remark 1. Defining a ratio 𝛽𝑖 via: 

𝛽𝑖 =
‖𝑣𝑖+1(𝑥)‖

‖𝑣𝑖(𝑥)‖
,                                                                                                                         (28) 

if this ratio stays less than one for all 𝑖, the approximate solution converges to the exact solution. 

 

5. Numerical results and discussions 

In this section, the methods presented in section three: the BOM, the shifted LOM, and the 

BrOM are applied to solve the nonlinear Blasius equation. However, first, the boundary 

conditions must be converted to initial conditions, i.e., the value of the second derivative at zero 

must be found. This was done in [27] using the homotopy-Padé method for approximate and 

finding that the best result is 𝑓′′(0) = 0.3320573. 

To solve this equation using the BOM method, we apply the technique from Section 3.1 with 

𝑛 = 3. Let us first assume the approximate solution as: 

 

 



IHJPAS. 37 (1) 2024 

365 
 

𝑓(𝑥) = 𝑐0 𝐵0,3(𝑥) + 𝑐1 𝐵1,3(𝑥) + 𝑐2 𝐵2,3(𝑥) + 𝑐3 𝐵3,3(𝑥) = 𝐶
𝑇𝜙(𝑥),                               (29) 

Where 𝐵0,3(𝑥) = 1 − 3𝑥 + 3𝑥
2 − 𝑥3, 𝐵1,3(𝑥) = 3𝑥 − 6𝑥

2 + 3𝑥3, 𝐵2,3(𝑥) = 3𝑥
2 − 3𝑥3,

𝐵3,3(𝑥) = 𝑥
3. 

From Eq.(8) we obtain the OM: 

𝐷𝐵 = [

−3 −1    0    0
   3 −1 −2    0
   0    2    1 −3
   0    0    1    3

] , (𝐷𝐵)
2 = [

   6   4    2    0
−12 −6    0    6
   6    0 −6 −12
   0    2    4    6

],  

(𝐷𝐵)
3 = [

−6 −6 −6 −6
   18    18    18    18
−18 −18 −18 −18
   6    6   6    6

]. 

After converting each 𝑓(𝑥) and its derivatives in the nonlinear Blasius equation and its initial 

conditions into terms of operational matrices, we substitute the collocation nodes Eq.(12) 

instead of each x, which gives us the following system of nonlinear algebraic equations: 

−6𝑐0 + 18𝑐1 − 18𝑐2 + 6𝑐3 +
1

2
(0.015625𝑐0 + 0.140625𝑐1 + 0.421875𝑐2 +

0.421875𝑐3)(1.5𝑐0 + 1.5𝑐1 − 7.5𝑐2 + 4.5𝑐3) = 0, 

−3𝑐0 + 3𝑐1 = 0,  

𝑐0 = 0, 

6𝑐0 − 12𝑐1 + 6𝑐2 = 0.3320573.                                                                                        (30) 

Finally, by solving (30) we obtain: 

𝑐0 = 0,   𝑐1 = 0,   𝑐2 = 0.055342883,   𝑐3 = 0.1635587513832345. 

So, the approximate solution will be as follows: 

𝑓(𝑥) = [0, 0, 0.055342883,      0.1635587513832345] [

1 − 3𝑥 + 3𝑥2 − 𝑥3

3𝑥 − 6𝑥2 + 3𝑥3

3𝑥2 − 3𝑥3

𝑥3

], 

= 0.16602865𝑥2 − 0.002469898616765498𝑥3. 

Thus, until 𝑛 = 11, the approximate solution will be: 

 

 𝑓(𝑥) = 0.16602865𝑥2 + 2.491375994395639 × 10−10𝑥3 − 1.878031063995422 ×

10−9𝑥4 − 0.000459416599124296𝑥5 − 2.537478849262697 × 10−8𝑥6 +

4.931977670707965 × 10−8𝑥7 + 0.000002433876389318357𝑥8 +

5.216138276864513 × 10−8𝑥9 − 2.538929066986384 × 10−8𝑥10 −

8.5669178417902 × 10−9𝑥11.                                                                                                                            (31) 

 



IHJPAS. 37 (1) 2024 

366 
 

Also, this equation is solved by using the shifted LOM method, and the description in Section 

3.2 is followed when 𝑛 = 3, the approximate solution is assumed to be as follows: 

𝑓(𝑥) = 𝑐0𝑃0(𝑥) + 𝑐1𝑃1(𝑥) + 𝑐2𝑃2(𝑥) + 𝑐3𝑃3(𝑥) = 𝐶
𝑇 𝜙(𝑥),                                            (32) 

Where 𝑃0(𝑥) = 1,   𝑃1(𝑥) = 2𝑥 − 1,   𝑃2(𝑥) = 6𝑥
2 − 6𝑥 + 1,  

𝑃3(𝑥) = 20𝑥
3 − 30𝑥2 + 12𝑥 − 1.  

From Eq.(19) we have the OM: 

𝐷𝐿 = [

0 0 0 0
2 0 0 0
0 6 0 0
2 0 10 0

] , (𝐷𝐿)
2 = [

0 0 0 0
0 0 0 0
12 0 0 0
0 60 0 0

] , (𝐷𝐿)
3 = [

0 0 0 0
0 0 0 0
0 0 0 0
120 0 0 0

]. 

After writing 𝑓(𝑥) and its derivatives in NODE and its initial conditions in terms of operational 

matrices, the collocation nodes Eq.(12) are substituted in place of each x to obtain the system 

of nonlinear algebraic equations: 

6𝑐0𝑐2 + 3. 𝑐1𝑐2 − 0.75𝑐2
2 + 120𝑐3 + 15. 𝑐0𝑐3 + 7.5𝑐1𝑐3 − 4.5𝑐2𝑐3 − 6.5625𝑐3

2 = 0, 

2𝑐1 − 6𝑐2 + 12𝑐3 = 0, 

𝑐0 − 𝑐1 + 𝑐2 − 𝑐3 = 0, 

12𝑐2 − 60𝑐3 = 0.3320573.                                                                                                  (33) 

Solving this system and substituting the value of 𝐶𝑇 into Eq.(32) yields the approximate 

solution: 

𝑓(𝑥) = −3.469446951953614 × 10−18 + 0.16602865𝑥2 − 0.002469898616765488𝑥3. 

 

Consequently, until 𝑛 = 11, we obtain the following approximate solution: 

𝑓(𝑥) = −3.696547211931118 × 10−12 + 3.253866343033706 × 10−10𝑥 +

0.16602864303933304𝑥2 + 6.283554605411232 × 10−8𝑥3 − 2.920702954934285 ×

10−7𝑥4 − 0.00045866085596424073𝑥5 − 0.000001148562469835698𝑥6 +

9.397408390670055 × 10−7𝑥7 + 0.000002163350548329646𝑥8.                              (34)  

 

Moreover, if we use the BrOM method, then the explanation in Section 3.3 is followed and the 

approximate solution is taken as assumed: 

𝑓(𝑥) = 𝑐0𝐵𝑟0(𝑥) + 𝑐1𝐵𝑟1(𝑥) + 𝑐2𝐵𝑟2(𝑥) + 𝑐3𝐵𝑟3(𝑥) = 𝐶
𝑇 𝐵𝑟(𝑥),                                (35) 

Where 𝐵𝑟0(𝑥) = 1,   𝐵𝑟1(𝑥) = 𝑥 −
1

2
,   𝐵𝑟2(𝑥) = 𝑥

2 − 𝑥 +
1

6
,   𝐵𝑟3(𝑥) = 𝑥

3 −
3

2
𝑥2 +

1

2
𝑥. 

Here, from Eq.(26) we have the OM: 

𝐷𝐵𝑟 = [

0 0 0 0
1 0 0 0
0 2 0 0
0 0 3 0

] , (𝐷𝐵𝑟)
2 = [

0 0 0 0
0 0 0 0
2 0 0 0
0 6 0 0

] , (𝐷𝐵𝑟)
3 = [

0 0 0 0
0 0 0 0
0 0 0 0
6 0 0 0

]. 



IHJPAS. 37 (1) 2024 

367 
 

After writing the nonlinear Blasius equation and its initial conditions in the form of operational 

matrices and putting the collocation nodes instead of each x, we obtain the following system: 

𝑐0𝑐2 + 0.25𝑐1𝑐2 − 0.02083333333333337𝑐2
2 + 6𝑐3 + 0.75𝑐0𝑐3 + 0.1875𝑐1𝑐3 −

0.06250000000000003𝑐2𝑐3 − 0.03515625𝑐3
2 = 0, 

𝑐1 − 𝑐2 +
𝑐3

2
= 0, 

𝑐0 −
𝑐1

2
+
𝑐2

6
= 0, 

2𝑐2 − 3𝑐3 = 0.3320573.                                                                                                      (36) 

Once this system is solved, the roots are substituted into Eq.(35), and we obtain the approximate 

solution as follows:  

𝑓(𝑥) = −6.93889390391 × 10−18 + 0.1660286500000003𝑥2 − 0.00246989861677𝑥3. 

 

Thus, up to 𝑛 = 11, the approximate solution will be: 

𝑓(𝑥) = −3.529427472531059 × 10−18 − 3.962758940434519 × 10−17𝑥 +

0.16602865000000003𝑥2 + 1.979266592913857 × 10−10𝑥3 − 1.499749921060817 ×

10−9𝑥4 − 0.000459418278828476𝑥5 − 2.065013702012278 × 10−8𝑥6 +

4.070826065187125 × 10−8𝑥7 + 0.000002444007699522198𝑥8 +

4.474523682223838 × 10−8𝑥9 − 2.232041556339086 × 10−8𝑥10 −

9.11508105682988 × 10−9𝑥11.                                                                                                                            

(37)  

The exact solution of the Blasius equation is unknown. Therefore, the maximum error 

remainder (𝑀𝐸𝑅𝑛) is calculated to determine the accuracy of the proposed methods. The 𝑀𝐸𝑅𝑛 

is given by: 

𝑀𝐸𝑅𝑛 = max
0≤𝑥≤1

|𝑓′′′(𝑥) +
1

2
𝑓(𝑥) 𝑓′′(𝑥)|.                                                                               (38) 

Figure 1 shows the logarithmic plots for 𝑀𝐸𝑅𝑛of the approximate solutions obtained by the 

three proposed methods at all 𝑛 iterations (𝑛 =  3 to 11). 

 

Figure 1.Logarithmic plots for the 𝑀𝐸𝑅𝑛 of BOM, LOM, and BrOM. 



IHJPAS. 37 (1) 2024 

368 
 

In addition, Table 1 shows the 𝑀𝐸𝑅𝑛 value for all 𝑛 studied in solving the nonlinear Blasius 

equation by using the proposed methods. 

 

Table 1: The comparison of 𝑀𝐸𝑅𝑛 between the BOM, LOM, and BrOM methods for 𝑛 = 3 to 11. 

𝑛 BOM LOM BrOM 

3 0.014819391700592988 0.014819391700592925 0.014819391700593073  

4 0.011427015738563973 0.011427015738558465 0.011427015738558968  

5 0.0004197059747635956 0.00041970597546251153 0.0004197059754639646  

6 0.00017127579193321196 0.00017127579194381904 0.00017127579189292255  

7 0.000025137988745926876 0.00002513798831121161 0.00002513798889495492  

8 0.000001263005001916894 0.000001263003365650816 0.000001263003561146565  

9 2.347967527072114

× 10−7 

5.772589885717771

× 10−7 

2.34800920268644 × 10−7  

10 1.588884090963915

× 10−8 

5.387294752107197

× 10−7 

1.59061837140361 × 10−8  

11 1.494825596637383

× 10−9 

5.385043960887126

× 10−7 

1.187559955162328 × 10−9  

 

From Figure 1 and Table 1, we can conclude that the value of the error decreases as 𝑛 increases. 

Thus, the BrOM method is better than the LOM method and slightly better than the BOM 

method. 

Moreover, a comparison was made between the approximate solutions found by the proposed 

methods at 𝑛 =  11, and the numerical solution obtained by the Range-Kutta method (RK4). 

This is shown in Figure 2, which shows a good agreement between them. 

            Figure 2. The comparison of the solutions of proposed methods and RK4 at 𝑛 = 11. 

 

To investigate the convergence of the solutions of the nonlinear Blasius equation, the 

convergence condition described in Section 4 is applied to the solutions of the proposed 

methods for all iterations n (n=3 to 11) by calculating the values of 𝛽𝑖 =
‖𝑣𝑖+1(𝑥)‖

‖𝑣𝑖(𝑥)‖
 , as shown in 



IHJPAS. 37 (1) 2024 

369 
 

Table 2, the values of 𝛽𝑖 for all 𝑖 ≥ 3 are less than 1, so these solutions converge to the exact 

solution. 

 

Table 2. The value of 𝛽𝑖 to the solutions of the proposed methods for 𝑛 = 3 to 11 when 𝑓′′(0) = 0.3320573. 

𝛽𝑖 BOM LOM BrOM 

𝛽3 0.016073750074326664 0.01607375007432271 0.016073750074321977 

𝛽4 0.29293061768966266 0.29293061766935835 0.2929306176693013 

𝛽5 0.03639714771755241 0.03639714778922806 0.03639714778673125 

𝛽6 0.24728635030717735 0.24728634929502963 0.24728634997162285 

𝛽7 0.09393295721721101 0.09393296065173386 0.09393296261024817 

𝛽8 0.0477950301697603 0.04779170199698017 0.047795085567504614 

𝛽9 0.13563608954043813 0.1356410881800995 0.1356460394028394 

𝛽10 0.047883452366514034 0.04874842641203115 0.04878982655985626 

On the other hand, in [5] the Blasius equation for 𝑓′′(0) = 1 was solved with BOM at 𝑛 = 11, 

so Figure 3 shows the difference between the approximate solution at  𝑓′′(0) = 0.3320573 

and 𝑓′′(0) = 1. 

Figure 3. Comparing the solution to the Blasius equation by BOM for two different values of 𝑓′′(0) at 𝑛 = 11. 

 

From the above figure, it can be seen that the solution is better at 𝑓′′(0) = 0.3320573. To 

illustrate this, the 𝑀𝐸𝑅𝑛 for both values of 𝑓′′(0) was calculated for all 𝑛 (𝑛 =  3 to 11) as in 

Table 3. 

 

 

 

 

 

 



IHJPAS. 37 (1) 2024 

370 
 

Table 3. The comparison of 𝑀𝐸𝑅𝑛 of the solutions of the Blasius equation by BOM at n=3 to 11 at two different 

values of 𝑓′′(0) 

𝑛 When 𝑓′′(0) = 0.3320573 When 𝑓′′(0) = 1 

3 0.014819391700592988 0.12364105477870169 

4 0.011427015738563973 0.09719389597650396 

5 0.0004197059747635956 0.011027011129465247 

6 0.00017127579193321196 0.004293535161145279 

7 0.000025137988745926876 0.0004907797800370872 

8 0.000001263005001916894 0.00009600784941188323 

9 2.347967527072114 × 10−7 0.00001612130638850573 

10 1.588884090963915 × 10−8 4.033942939685175 × 10−7 

11 1.494825596637383 × 10−9 2.636242850684311 × 10−7 

 

Thus, it can be clearly seen that 𝑀𝐸𝑅𝑛 converges to zero faster in the case of 𝑓′′(0) =

0.3320573 than in the case of at 𝑓′′(0) = 1, as 𝑛 increases. 

Moreover, we will solve this problem with the initial condition 𝑓′′(0) = 1 using the shifted 

LOM and the BrOM. The solution with the shifted LOM when 𝑛 = 11 is as follows: 

𝑓(𝑥) = 2.717523012968703 × 10−11 − 2.989454184199047 × 10−9𝑥 +

0.5000000807212571𝑥2 − 8.890725537863409 × 10−7𝑥3 +

0.00000538235410165247𝑥4 − 0.00418570492658484𝑥5 +

0.000041144785309415706𝑥6 − 0.000054417856055302455𝑥7 +

0.00010990985811043829𝑥8 − 0.000015116787731315758𝑥9.                                   (39) 

If the BrOM apply, the approximate solution is as follows: 

𝑓(𝑥) = 4.28022688906543 × 10−17 + 4.126066892665148 × 10−17𝑥 + 0.5𝑥2 +

4.379568419265968 × 10−8𝑥3 − 3.321493629909657 × 10−7𝑥4 −

0.0041651294231496725𝑥5 − 0.000004588795148760594𝑥6 +

0.000009070467696073757𝑥7 + 0.000056311613028286645𝑥8 +

0.000010056188911186578𝑥9 − 0.000005051467149564113𝑥10 +

5.856448890320508 × 10−9𝑥11.                                                                                        (40) 

Table 4 shows a comparison between the 𝑀𝐸𝑅𝑛 of the three proposed methods for this problem 

for all iterations n (n=3 to 11). It can be seen that the accuracy increases as n increases. Also, 

the BrOM and the BOM methods are better than the LOM method, and the BrOM method is 

slightly better than the BOM method. 

 

 

 

 



IHJPAS. 37 (1) 2024 

371 
 

Table 4. The comparison of 𝑀𝐸𝑅𝑛 between the BOM, LOM, and BrOM methods for 𝑛 = 3 to 11 at 𝑓′′(0) = 1. 

𝑛 BOM LOM BrOM 

3 0.12364105477870169 0.12364105477870149 0.12364105477870134 

4 0.09719389597650396 0.09719389597647592 0.09719389597650352 

5 0.011027011129465247 0.011027011129505574 0.011027011129506844 

6 0.004293535161145279 0.004293535161116567 0.00429353516106645 

7 0.0004907797800370872 0.0004907797800190347 0.0004907797804556991 

8 0.00009600784941188323 0.00009600784484897892 0.00009600784504393617 

9 0.00001612130638850573 0.000016121319056666043 0.000016121340651173055 

1

0 

4.033942939685175

× 10−7 

0.000005672419462032785 4.034414582290297

× 10−7 

1

1 

2.636242850684311

× 10−7 

0.000005624000721099476 2.627741051773592

× 10−7 

Moreover, the convergence condition given in Section 4 is applied to the solutions of the 

proposed methods for all iterations n (n=3 to 11) when 𝑓′′(0) = 1 by calculating the values of 

𝛽𝑖 =
‖𝑣𝑖+1(𝑥)‖

‖𝑣𝑖(𝑥)‖
 , as shown in Table 5, the values of 𝛽𝑖 for all 𝑖 ≥ 3 are less than 1. Consequently, 

the convergence condition is satisfied. 

 

Table 5. The value of 𝛽𝑖 to the solutions of the proposed methods for 𝑛 = 3 to 11 when 𝑓′′(0) = 1. 

𝛽𝑖 BOM LOM BrOM 

𝛽3 0.04510659083967248 0.04510659083966797 0.04510659083967211 

𝛽4 0.277633864539339 0.27763386453933614 0.04510659083967211 

𝛽5 0.11985758340437944 0.11985758340397297 0.1198575834036301 

𝛽6 0.23281227458321288 0.23281227457927967 0.23281227459976303 

𝛽7 0.06500373117417597 0.06500373174199665 0.0650037318024709 

𝛽8 0.21189027133781216 0.21189028400881857 0.21189031026366478 

𝛽9 0.1197929121542225 0.11979856224759194 0.11979340082544765 

𝛽10 0.008393018377721006 0.008431707509080347 0.008431055448947825 

 

6. Conclusions 

In this work, the operational matrices of differentiation were used based on: Bernstein, shifted 

Legendre, and Bernoulli polynomials to solve the nonlinear Blasius equation. We concluded 

that the Bernoulli operational matrix method is better than the Bernstein operational matrix 

method and the shifted Legendre operational matrix method.  We also compared the results 

numerically with the Runge-Kutta method and found good agreement between them. The 

calculations in this study were performed using the Mathematica® 12 program. Also, the 

maximal error remainder value for the proposed approximation methods was calculated. 



IHJPAS. 37 (1) 2024 

372 
 

Moreover, after examining a certain value of the second derivative (𝑓′′(0) = 0.3320573) 

calculated in [27], the result was compared with the value obtained in [5] (𝑓′′(0) = 1). It was 

found that the solution is more accurate at 𝑓′′(0) = 0.3320573. 

Acknowledgment 

The authors are greatly appreciated the referees for their valuable comments and suggestions 

for improving the paper 

Conflict of Interest  

The authors declare that they have no conflicts of interest. 

Funding 

There is no financial support in preparation for the publication. 

References 

1. Murphy, G.M. Ordinary Differential Equations and Their Solutions; Dover Publications, 

Inc., New York, 1960. 

2. Boyce, W.E.; DiPrima, R.C. Elementary differential equations and boundary value 

problems, Ed.;9th ed. 2009; John Wiley & Sons, Inc., United States of America, 2009; ISBN 

9780470383346. 

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

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

https://doi.org/10.1080/00207720903154783 

4. Pirabaharan, P.; Chandrakumar, R.D.; A computational method for solving a class of 

singular boundary value problems arising in science and engineering, Egyptian Journal of 

Basic and Applied Sciences 2016, 3(4), 383-391. DOI: 

https://doi.org/10.1016/j.ejbas.2016.09.004 

5. Khataybeh, S.; 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. DOI: 

https://doi.org/10.1016/j.jksus.2018.05.002 

6. Asgari, M.; Ezzati, R.; Using operational matrix of two-dimensional Bernstein polynomials 

for solving two-dimensional integral equations of fractional order, Applied Mathematics 

and Computation 2017, 307(C), 290-298. DOI: https://doi.org/10.1016/j.amc.2017.03.012 

7. Hesameddini, E.; Shahbazi, M.; Two-dimensional shifted Legendre polynomials 

operational matrix method for solving the two- dimensional integral equations of fractional 

order, Applied Mathematics and Computation 2018, 322(C), 40-54. DOI: 

https://doi.org/10.1016/j.amc.2017.11.024 

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

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

https://doi.org/10.1515/nleng-2017-0157 

9. Bazm, S.; Bernoulli polynomials for the numerical solution of some classes of linear and 

nonlinear integral equations, Journal of Computational and Applied Mathematics 2015, 

275(C), 44-60. DOI: https://doi.org/10.1016/j.cam.2014.07.018 

https://doi.org/10.1080/00207720903154783
https://doi.org/10.1016/j.ejbas.2016.09.004
https://doi.org/10.1016/j.jksus.2018.05.002
https://doi.org/10.1016/j.amc.2017.03.012
https://doi.org/10.1016/j.amc.2017.11.024
https://doi.org/10.1515/nleng-2017-0157
https://doi.org/10.1016/j.cam.2014.07.018


IHJPAS. 37 (1) 2024 

373 
 

10. Al-Jawary, M.A.; Ibraheem, G.H.; Tow meshless methods for solving nonlinear ordinary 

differential equations in engineering and applied sciences, Nonlinear Engineering 2020, 

9(1), 244-255. DOI: https://doi.org/10.1515/nleng-2020-0012 

11. 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: https://doi.org/10.1016/j.aej.2020.07.008 

12. 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. DOI: https://doi.org/10.1080/16583655.2019.1580662 

13. Bani-Ahmad, F.; Alomari, A.K.; Bataineh, A.S.; Sulaiman, J.; Hashim, I.; On the 

approximate solutions of systems of ODEs by Legendre operational matrix of 

differentiation, Italian Journal of Pure and Applied Mathematics 2016, 36, 483-

494.https://ijpam.uniud.it/online_issue/201636/42 

14. Kumar, S.; Pandey, P.; Das, S.; Craciun, E.-M.; Numerical solution of two dimensional 

reaction-diffusion equation using operational matrix method based on Genocchi 

polynomial-Part I: Genocchi polynomial and opperatorial matrix, Proceedings of the 

Romanian Academy, Series A 2019, 20(4), 393-399. 

https://acad.ro/sectii2002/proceedings/doc2019-4 

15. Loh, J.R.; Phang, C.; Numerical solution of Fredholm fractional integro-differential 

equation with right-sided Caputo’s derivative using Bernoulli polynomials operational 

matrix of fractional derivative, Mediterranean Journal of Mathematics 2019, 16(2), 1-25. 

 DOI: https://doi.org/10.1007/s00009-019-1300-7 

16. Zeghdane, R.; Numerical solution of stochastic integral equations by using Bernoulli 

operational matrix, Mathematics and Computers in Simulation 2019, 165(C), 238-254. 

DOI: https://doi.org/10.1016/j.matcom.2019.03.005 

17. Jasim, S.M.N.; Ibraheem, G.H.; Fractional Pantograph Delay Equations Solving by the 

Meshless Methods, Ibn AL-Haitham Journal For Pure and Applied Sciences 2023, 36, 382-

397. DOI: https://doi.org/10.30526/36.3.3076 

18. Alshbool, M.H.T.; Mohammad, M.; Isik, O.; Hashim, I.; Fractional Bernstein operational 

matrices for solving integro-differential equations involved by Caputo fractional derivative, 

Results in Applied Mathematics 2022, 14, 100258. DOI: 

https://doi.org/10.1016/j.rinam.2022.100258 

19. Salih, A.A.; Shihab, S.; New operational matrices approach for optimal control based on 

modified Chebyshev polynomials, Samarra Journal of Pure and Applied Science 2020, 

2(2), 68–78. DOI: http://dx.doi.org/10.54153/sjpas.2020.v2i2.115 

20. Jalal, R.; Shihab, S.; Abed Alhadi, M.; Rasheed, M.; Spectral Numerical Algorithm for 

Solving Optimal Control Using Boubaker-Turki Operational Matrices, Journal of Physics: 

Conference Series, IOP Publishing 2020, 1660(1), 012090. DOI: 

https://doi.org/10.1088/1742-6596/1660/1/012090 

21. 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. DOI: https://doi.org/10.1016/j.rinp.2019.102836 

22. Kaur, H.; Mishra, V.; Mittal, R.C.; Numerical solution of a laminar viscous flow boundary 

layer equation using uniform Haar wavelet quasi-linearization method, International 

https://doi.org/10.1515/nleng-2020-0012
https://doi.org/10.1016/j.aej.2020.07.008
https://doi.org/10.1080/16583655.2019.1580662
https://ijpam.uniud.it/online_issue/201636/42-Bani-AhmadAlomariBatainehSulaimanHashim.pdf
https://acad.ro/sectii2002/proceedings/doc2019-4/10-Craciun.pdf
https://doi.org/10.1007/s00009-019-1300-7
https://doi.org/10.1016/j.matcom.2019.03.005
https://doi.org/10.30526/36.3.3076
http://dx.doi.org/10.1016/j.rinam.2022.100258
http://dx.doi.org/10.54153/sjpas.2020.v2i2.115
https://ui.adsabs.harvard.edu/link_gateway/2020JPhCS1660a2090J/doi:10.1088/1742-6596/1660/1/012090
https://doi.org/10.1016/j.rinp.2019.102836


IHJPAS. 37 (1) 2024 

374 
 

Journal of Mathematical and Computational Sciences 2013, 79, 1410-1415. DOI: 

http://dx.doi.org/10.5281/zenodo.1087368 

23. Cortell, R.; Numerical solutions of the classical Blasius flat-plate problem, Applied 

Mathematics and Computation 2005, 170(1), 706-710. DOI: 

https://doi.org/10.1016/j.amc.2004.12.037 

24. [24] He, J.; Approximate analytical solution of Blasius’ equation, Communications in 

Nonlinear Science and Numerical Simulation 1999, 4(1), 75-78. DOI: 

https://doi.org/10.1016/S1007-5704(99)90063-1 

25. 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. DOI: https://doi.org/10.1016/j.chaos.2005.10.071 

26. Aminikhah, H.; An analytical approximation for solving nonlinear Blasius equation by 

NHPM, Numerical Methods for Partial Differential Equations 2010, 26(6), 1291-1299. 

DOI: https://doi.org/10.1002/num.20490 

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

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

2003-2016. DOI: https://doi.org/10.1016/j.cnsns.2009.09.002 

28. Mohammed Ali, M.N.; A New operational matrix of derivative for orthonormal Bernstein 

polynomial’s, Baghdad Science Journal 2014, 11(3), 1295-1300. 

DOI:https://doi.org/10.21123/bsj.2014.11.3.1295-1300  

29. 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-1281. 

DOI: https://doi.org/10.21123/bsj.2014.11.3.1274-1281  

30. Turkyilmazoglu, M.; Convergent optimal variational iteration method and applications to 

heat and fluid flow problems, International Journal of Numerical Methods for Heat & Fluid 

Flow 2016, 26(3/4), 790–804. DOI: https://doi.org/10.1108/HFF-09-2015-0353 

 

 

 

http://dx.doi.org/10.5281/zenodo.1087368
https://doi.org/10.1016/j.amc.2004.12.037
https://doi.org/10.1016/S1007-5704(99)90063-1
https://doi.org/10.1016/j.chaos.2005.10.071
https://doi.org/10.1002/num.20490
https://doi.org/10.1016/j.cnsns.2009.09.002
https://doi.org/10.21123/bsj.2014.11.3.1295-1300
https://doi.org/10.21123/bsj.2014.11.3.1274-1281
https://doi.org/10.1108/HFF-09-2015-0353