388 

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

IHJPAS. 37 (2) 2024  

Ibn Al-Haitham Journal for Pure and Applied Sciences 

Journal homepage: jih.uobaghdad.edu.iq 

PISSN: 1609-4042, EISSN: 2521-3407 

 

 

 

 

Firas A. Fawzi*1  , Hadeer M. Globe2  and Nizam G. Ghawadri 3  
1,2 Department of Mathematics, Faculty of Computer Science and Mathematics, Tikrit University,  

Salah Al-Din, Iraq. 
3 Doctor of Mathematics, Ministry of Education, Jenin, Palestine. 

*Corresponding Author. 

 

 

 

Abstract   
 

This paper presents two important contributions to the field of numerical analysis for third-order 

ordinary differential equations (ODEs). First, a new class of direct implicit Runge-Kutta (RK) 

processes, called RKTDIO, is introduced as solutions to third-order ODEs. Secondly, it develops 

the ERKTDIO method, which is an embedded pairwise diagonal implicit RK method. The study 

begins by introducing the theory of relevant-colored trees and B-series as fundamental concepts. 

By utilizing the order constraints, two RKTDIO methods are derived: a fifth-order method with 

three stages and a sixth-order method with four stages. In addition, an embedded method called 

ERKTDIO6(5) is derived, which has orders six and five. The derivation of the embedded method 

includes strategies to ensure that the higher-order method achieves high accuracy while the lower-

order method provides optimal error estimates. To evaluate the effectiveness of the proposed 

methods, variable step-size codes are developed and applied to a set of specific third-order 

problems. The numerical evaluation involves converting the problems into a system of first-order 

ODEs and comparing the results with existing methods in terms of accuracy and function 

evaluations. The numerical demonstrations emphasise the superior performance and efficiency of 

the new methods in solving third-order ODEs. The comparative analysis shows the accuracy 

achieved by the higher-order method and the improved error estimation of the lower-order method. 

The results validate the efficacy of the proposed approaches and their potential for practical 

applications in various domains.  

Keywords:  Third-order ODEs, order conditions, B-series, Relevant-colored trees.  
 

1. Introduction 

Third-order differential equations can be found in various fields, such as applied sciences, neural 

network engineering [1,2], fluid dynamics [3] and thin film flow [4]. The aim of this paper is to 

develop and explain a computational method for solving initial value problems of third-order 

differential equations. 

𝛼′′′(𝑥)  =  𝜇 (𝑥, 𝛼(𝑥)),          𝑥 ≥  𝑥0,                                                                                           (1) 

With initial conditions 

α(xn) = αn  ,   β(xn) = 𝛼𝑛
′  ,    γ(xn) = α𝑛

′′ . 

Received: 13 April 2023               Accepted: 12 June 2023                      Published: 20 April 2024 

 

 

Efficient Embedded Diagonal Implicit Runge-Kutta Method for Directly 

Solving Third Order ODEs 

 

doi.org/10.30526/37.2.3409 

 

https://creativecommons.org/licenses/by/4.0/
https://jih.uobaghdad.edu.iq/index.php/j/index#1609-4042
https://jih.uobaghdad.edu.iq/index.php/j/index#2521-3407
https://orcid.org/my-orcid?orcid=0000-0003-0939-7940
mailto:firasadil01@tu.edu.iq
https://orcid.org/0009-0008-2163-2484
mailto:hadeermohamadjalob@gmail.com
https://orcid.org/0000-0002-4335-1789
mailto:nizamghawadri@gmail.com


IHJPAS. 37 (2) 2024 

 

389 
 

 

Where (𝑥)∈ 𝑅𝑑,  𝑓:  𝑅𝑑× 𝑅𝑑 → 𝑅𝑑 has a continuous value but lacks the first and second derivatives. 

Numerical solutions are often necessary for third-order differential equations since analytical 

solutions are often not available. Some researchers have used classical methods to solve higher-

order differential equations by converting them into a system of first-order differential equations 

[5-16,27]. However, this method can be computationally intensive and time consuming. Direct 

numerical approaches have been proposed to reduce the computation time, but they require a 

second method to obtain initial values for the numerical solutions. Implicit methods are useful 

because they can achieve high accuracy with fewer steps, making it easier to find solutions to 

complex problems. Several researchers have developed embedded Runge-Kutta methods with 

high algebraic orders for solving third-order differential equations [17-18]. Ismail et. al [19] 

proposed the Singly Embedded Diagonally Implicit Runge-Kutta (SDIRK) method to combine 

delay differential equations (DDEs) and compared the computational results. The researchers in 

[20,21,26] developed a novel embedded explicit and implicit Runge-Kutta method for solving 

special third and fourth-order problems. 

The main objective of this work is to develop a new method called RKTDIO, an implicit one-step 

Runge-Kutta method designed for directly solving specific third-order differential equations. This 

method is developed using the theory of relevant-colored trees theory and also involves the 

derivation of embedded diagonally implicit Runge-Kutta methods for the direct integration of 

certain third-order differential equations. 

The structure of this paper is as follows: Section 2 provides the formulation concept for the 

RKTDIO approach for the direct integration of certain third-order ODEs. In Section 3, we develop 

the novel theory of relevant-colored trees theory and the corresponding B-series theory. Section 4 

contains the derivation of the ordering criteria of the RKTDIO method. Section 5 presents the 

construction of a three-stage RKTDIO approach for order five and a four-stage method for order 

six. Section 6 describes the derivation of the embedded ERKTDIO6(5) method. To demonstrate 

the efficiency and effectiveness of the RKTDIO and ERKTDIO6(5) methods compared to the 

methods currently used in the scientific literature, we give numerical findings in Section 7. Section 

8 provides conclusions. 

 

2.  The formulation of RKTDIO method 

By turning it into a system of first-order ODEs, the special third-order IVP (1) can be solved as 

follows: 

(

α(x)

β(x)

γ(x)
)

′

= (

β(x)

γ(x)

μ(x, α(x))
)                                                                                                                         (2) 

With initial conditions 

α(xn) = αn  ,   β(xn) = 𝛼𝑛
′  ,    γ(xn) = α𝑛

′′. 

The Runge-Kutta method first-order is used to obtain the following system of equations 

αi= αn + h ∑ aijα′j
s
j=1 ,                                                                                                                               (3)                                                                                               

α′i = 𝛼𝑛
′ + h ∑ aijα′′j ,

s
j=1                                                                                                                          (4)                                                                                             

 α𝑖
′′ = α𝑛

′′ + h ∑ aij μ(xn + cjh, αj) ,s
j=1                                                                                                  (5)                                                           

αn+1 =  αn + h ∑ biα𝑖
′    ,                                                                                                                      (6)s

i=1        



IHJPAS. 37 (2) 2024 

 

390 
 

α𝑛+1
′ =  α𝑛

′ + h ∑ biα𝑖
′′s

i=1   ,                                                                                                                     (7)                                                                                    

α𝑛+1
′′ =  α𝑛

′′ + h ∑ biμ
s
i=1 (xn + cih, αi) .                                                                                               (8)                                                                

 

When we ignore  𝛼𝑖
′, 𝛼𝑖

′′ and α𝑖
′′′ , from (3) – (8) we conclude 

αi = αn + h ∑ aij

s

j=1

α𝑛
′ + h2 ∑ aij

s

j,k=1

αjkα𝑛
′′

+  h3 ∑ aijajkakl μ(xn + cih, αl), i = 1, … , s.                                                   (9)

s

j,k,l

 

𝛼𝑛+1 = 𝛼𝑛 + ℎ ∑ 𝑏𝑖𝛼𝑛
′ + ℎ2 ∑ 𝑏𝑖

𝑠

𝑖,𝑗=1

𝑎𝑖𝑗𝛼𝑛
′′

𝑠

𝑖=1

+ ℎ3 ∑ 𝑏𝑖𝑎𝑖𝑗𝑎𝑗𝑘  𝜇(𝑥𝑛 + 𝑐𝑘ℎ,

𝑠

𝑖,𝑗,𝑘=1

𝛼𝑘),                                                                       (10) 

𝛼′𝑛+1 = 𝛼𝑛
′ + ℎ ∑ 𝑏𝑖𝛼𝑛

′′ + ℎ2 ∑ 𝑏𝑖𝑎𝑖𝑗 𝜇𝑠
𝑖.𝑗=1 (𝑥𝑛 + 𝑐𝑗ℎ𝑠

𝑖=1 , 𝛼𝑗),                                                       (11)  

𝛼𝑛+1
′′ =  𝛼𝑛

′′ + h ∑ bi μ(xn + cih,

s

i=1

 αi)  .                                                                                              (12) 

We assume that 

∑ aij = ci,

s

j=1

        ∑ aijajk =
1

2

s

j,k=1

 ci
2 ,         ∑ bi = 1

s

i=1

,          ∑ biaij =
1

2

s

i,j=1

, 

∑ biaij =  b𝑖
′′  ,     ∑ bjajkakl = b𝑖

′′

s

j,k=1

s

j=1

, ∑ aik

s

k,l,r=1

aklalr = �̂�𝑖𝑗 , 𝑖 = 1, … , 𝑠. 

As a result, we are able to solve the specific third-order IVP (1), indicated by the RKTDIO 

approach, using the following direct integration method. Thus, the following formula is used to 

represent the s-stage RKTDIO approach for the numerical solution of the IVP (1): 

αn+1 =  αn + h𝛼𝑛
′ +  

1

2
 h2𝛼𝑛

′′ + h3  ∑ bi

s

i=1

 μ(xn + cih, αi),                                                            (13) 

𝛼′𝑛+1 =  𝛼𝑛
′ + h𝛼𝑛

′′ + h2 ∑ b𝑖
′ μ(xn + cih,

s

i=1

 α𝑖
′) ,                                                                              (14) 

𝛼𝑛+1
′′ = 𝛼𝑛

′′ + h ∑ b𝑖
′′

s

i=1

 μ (xn + cih, α𝑖
′) ,                                                                                             (15)  

αi =  αn + cih𝛼𝑛
′ +  

1

2
ci

2h2𝛼𝑛
′′ + h3 + ∑ �̂�ij   μ(xn + cih,

s

i,j=1

 α𝑗
′).                                                  (16) 

All RKTDIO parameters  aij, bi, b𝑖
′, b𝑖

′′ , and  ci are real numbers and i, j = 1, 2, . . . , s. 

The RKTDIO method (13)–(16) can be expressed in Butcher tableau as follows 

 

 



IHJPAS. 37 (2) 2024 

 

391 
 

 

 

 

 

 

: 

 

3. B-series and associated relevant-colored trees 

The essential definitions of lemmas and associated theorems that are utilized throughout this article 

will be covered in this part. 

Definition 3.1 [22,23]: Order of the RKTDIO Method for Third-Order ODEs 

The RKTDIO method (13)-(16) has order 𝑝 when the third-order ODE (1) with the assumption 

α(xn) = αn,   α′(xn) = α′
n, α′′(xn) =  α′′

n   hence the local truncation error norms of the exact 

solution and the first, and second derivatives of the solution must satisfied. 

α(xn + h) − αn+1 = O(hp+1), 

α′(xn + h) − α′n+1 = O(hp+1), 

α′′(xn + h) − α′′
n+1 = O(hp+1),                                                                                                        (17) 

The following autonomous form of third-order IVP must be used in order to derive the algebraic 

order conditions for the RKTDIO technique (13)–(16). 

α(3)(x) = ρ(α(x)),                                                                                                                                 (18) 

With initial conditions 

α(xn) = αn,    α′(xn) = α′
n,     α′′(xn) = α𝑛

′′. 
By extending IVP (1) with a one-dimensional vector 𝑤 =  𝑥, the autonomous problem can be 

expressed equivalently to the third-order initial value problem (1) as follows: 

w3 = 0,                                                                                                                                                      (19) 

α(3) = ρ(w, α),                                                                                                                                         (20) 

w(xn) = wn    w′(xn) = w′
n = 1,      w′′(xn) = w′′

n = 0,                                                          (21)   

α(xn) = αn,        α′(xn) = α′
n,         α′′(xn) = α′′n.                                                                         (22) 

Applying RKTDIO method (13)–(16) to the scheme (19)–(22), we obtain 

Wi = wn + cih w′
n +

1

2
ci

2h2 w′′
n ,                                                                                                   (23)  

αi = αn + cihα′
n +

1

2
ci

2h2 α′′
n + h3 ∑ aij μ(Wj, αj)  ,

       

s

i,j=1

                                                              (24) 

wn+1 = wn + h w′
n +

1

2
h2w′′

n  ,                                                                                                       (25) 

w′
n+1 = w′

n + h w′′
n  ,                                                                                                                        (26)    

w′′
n+1 = w′′

n      ,                                                                                                                                   (27)  

αn+1 = αn + h α′
n +

1

2
h2α′′

n + h3 ∑ bi μ(Wi

s

i=1

, αi),                                                                     (28)   

α′n+1 = α′
n + h α′′

n + h2 ∑ b′
i μ(Wi,

s

i=1

 αi)   ,                                                                                  (29)   

𝑐1

⋮
𝑐𝑠

 
�̂�11 … �̂�1𝑠

⋮ ⋱ ⋮
�̂�𝑠1 … �̂�𝑠𝑠

 

 𝑏1 … 𝑏𝑠 

𝑏1
′ … 𝑏𝑠

′  

𝑏1
′′ … 𝑏𝑠

′′ 



IHJPAS. 37 (2) 2024 

 

392 
 

α′′
n+1 =  α′′

n + h ∑ b′′
i 

s

i=1

μ(Wi, αi).                                                                                                  (30) 

Substituting Eq. (21) into system of Eqs. (23)–(30), we get 

Wi = xn + cih,                                                                                                                                          (31)        

wn+1 = xn + h,                                                                                                                                        (32)                                                                                                   

w′
n+1 = 1,                                                                                                                                                (33)                                                                                                       

w′′
n+1 = 0,                                                                                                                                               (34)                                                                                                                                         

αn+1 = αn + hα′
n +

1

2
 h2α′′

n + h3 ∑ bi
s
i=1  μ(xn + cih, αi),                                                        (35)            

α′
n+1 = α′

n + hα′′
n + h2 ∑ b′

i
s
i=1  μ(xn + cih, αi),                                                                (36) 

α′′
n+1 = α′′

n + h ∑ b′′
i

s
i=1  μ(xn + cih, αi),                                                                             (37) 

αi = αn + cihα′
n +

1

2
ci

2h2α′′
n + h3 ∑ �̂�𝑖𝑗 μ(xn + cih

s
i,j=1 , αj).                                             (38) 

We find that the system of equations (13)–(16), which is generated by using the RKTDIO approach 

on the non-autonomous problem (1), is totally similar to Eqs. (35)– (38). Hence, discussing the 

numerical solutions of autonomous form is sufficient (18). Hence, the RKTDIO method (13)–(16) 

can be implemented as follows. 

αn+1 = αn + h α′
n +

1

2
 h2α′′

n + h3 ∑ bi

s

i=1

 μ(αi), 

α′
n+1 = α′

n + hα′′
n + h2 ∑ b′

i

s

i=1

 μ(αi), 

α′′
n+1 = α′′

n + h ∑ b′′
i

s

i=1

 μ(αi), 

αi = αn + cihα′
n +

1

2
 ci

2h2α′′
n + h3 ∑ �̂�𝑖𝑗

s
i,j=1 μ(αj).                                                               (39) 

The following elementary differentials are obtained by using the elementary differential notation 

on the analytical solution 𝛼(𝑥). [15] 

α(1) = α′,            α(2) = α′′,        α(3) = μ,       

α(4) = μ′α′,         α(5) = μ′′(α′, α′) + μ′α′′(𝛼′, 𝛼′),        

α(6) = 3𝛼′′(𝛼′, 𝛼′′) + μ′′′(𝛼′, 𝛼′, 𝛼′) + μ′α′′.                                                                            (40) 

These processes very quickly get more difficult as the order increases. The optimum method for 

overcoming this challenge, according to [25], will be to use a graphical representation with a few 

modifications for third-order ODEs, denoted by relevant-colored trees. The three sorts of nodes in 

the relevant-colored trees are "meager," "black ball", and "white ball", and they are connected by 

arcs. In these trees, we specifically use the end meager node to denote each 𝛼′, the end black ball 

node to denote each 𝛼′′, the end white ball to denote each 𝜇, and each arc to denote each arc, 

leaving this node to represent the 𝑚 − 𝑡ℎ derivative of 𝜇 with respect to 𝛼. In addition, the symbols 

𝑡1 and 𝑡2 denote the first-order and the second-order tree and 𝑡3 the third-order tree, respectively 

(see Figure 1) 

 



IHJPAS. 37 (2) 2024 

 

393 
 

𝜏1 =  
𝜏2 =  

𝜏3 =  

 

Figure 1. The relevant-colored trees 

 

Here, we'll go through some essential definitions for the relevant-colored trees and related B-series 

that are necessary for this work. 

Definition 3.2 [22,23]: Relevant-Colored Trees (RT) and Meager Node 

The following definitions are repeated for the set of relevant-colored trees (RT): 

1- The trees 𝑡1, 𝑡2 and  𝑡3 above are all in RT, and the tree 𝑡1 includes only one meager node 

(known as the root). 

2- If 𝜏1, 𝜏2, . . . , 𝜏𝑚  ∈ 𝑅𝑇, then 𝜏 =  [𝜏1, 𝜏2, … , 𝜏𝑚]3 is the tree obtained by linking the roots 

𝜏1, 𝜏2, . . . , 𝜏𝑚  and the root is the ‘‘meager node’’ 𝑡1 is at the bottom. The subscript 3 is to 

mention that the trees of the roots of 𝜏1, 𝜏2, . . . , 𝜏𝑚 onto the tree 𝑡3 contain a chain of three 

nodes. 

We will employ the following principles to create the appropriate colored trees: 

1- The ‘‘meager’’ node is always the root. 

2- A ‘‘meager’’ node has only single child and that child should be ‘‘black ball’’. 

3- A ‘‘black ball’’ node has only single child and that child should be ‘‘white ball’’. 

Definition 3.3 [13,14]: Order Function 𝝆(𝝉) for Relevant-Colored Trees (RT) 

The order 𝜌(𝜏) and symmetry 𝜎(𝜏) functions are defined recursively as follows: 

1-  𝜌(𝑡1) =  1, 𝜌(𝑡2) =  2, 𝜌(𝑡3)  =  3, 

2- 𝜎(𝑡1)  =  1, 𝜎(𝑡2)  =  1, 𝜎(𝑡3)  =  1,  

3- If 𝜏 =  [𝜏1, 𝜏2, … , 𝜏𝑚]3 for each 𝜏 ∈  𝑅𝑇, then ρ(τ ) = 3 + ∑ ρ(τi)
m
i=1  and 𝜎(𝜏 )  =

∏ 𝜎(𝜏𝑖)(𝑣1!  𝑣2!𝑚
𝑖=1 … ) where 𝜌(𝜏) is the number of nodes of 𝜏, ∀τ ∈  RT and  ν1!  ν2! …  

count equal trees among 𝜏1, 𝜏2, . . . , 𝜏𝑚.  Then we can define the set 𝑆𝑟 which consist of 

every trees 𝑅𝑇 of order 𝑟. 

Definition 3.4 [22,23]: Elementary Differential and B-Series on Relevant-Colored Trees (RT) 

for the RKTDIO Approach 

The elementary differential for every tree 𝜏 ∈  𝑅𝑇 is a function F(τ): Rd × Rd × Rd → Rd , 

recursively defined on RT as follows 

1- 𝕌(𝑡1)(𝛼, 𝛼′, 𝛼′′) = 𝛼′, 𝕌(𝑡2)(𝛼, 𝛼′, 𝛼′′) = 𝛼′′, 𝕌(𝑡3) (𝛼, 𝛼′, 𝛼′′) = 𝜇(𝛼), 

2- 𝕌 (𝜏 )  =  𝜇(𝑚)(𝛼) (𝕌 (𝜏1)(𝛼, 𝛼′, 𝛼′′), . . . , 𝕌 (𝜏𝑚)(𝛼, 𝛼′, 𝛼′′)) for 𝜏 =  [𝜏1, 𝜏2, … , 𝜏𝑚]3. 

  We expanded these definitions to provide the definition of B-series on the set RT of the relevant-

colored trees for the RKTDIO approach, which was motivated by the definitions of B-series on 

the root trees in [23] and the tri-colored trees in [24]. 

Definition 3.5 [22,23]: B-Series Representation 

Let δ: RT ∪ {∅} → Rd be a mapping, then we can give the form of a formal series as follows: 

B(δ, y) = δ(∅) + ∑
hρ(τ)

σ(τ)
δ(τ)τ∈RT  𝕌(τ)(α, α′α′′) ,                                                                   (41) 

Which is called the B-series. 

 



IHJPAS. 37 (2) 2024 

 

394 
 

      We present the following crucial lemma that is important to this derivation in order to 

accomplish the main goal of this research, which is the derivation of the order conditions of the 

RKTDIO technique. 

Lemma 3.1 [22,23]: Let 𝛿 be a function 𝛿 ∶ 𝑅𝑇 ∪ {∅} → 𝑅𝑑 with 𝛿(∅) = 1. Thus ℎ3𝜇(𝐵(𝛿, 𝛼)) 

is also a B-series ℎ3𝜇(𝐵(𝛿, 𝛼)) = 𝐵(𝛿′, 𝛼) where 𝛿′(∅) = 0,   𝛿′(𝑡1) = 0,   𝛿′(𝑡2) =

0,   𝛿′(𝑡3) = 1, and for 𝜏 =  [𝜏1, 𝜏2, … , 𝜏𝑚]3 ∈ 𝑅𝑇 ,   𝛿′(𝜏) = 𝛿(𝜏1 ),.  .  . 𝛿(𝜏𝑚). 

Lemma 3.2 [22,23]: If we suppose that the analytic solution of (18) is a B-series 𝐵(𝜗, 𝛼0) with a 

real function  𝜗 which is defined on 𝑅𝑇 ∪ {∅}, then 

 𝜗(∅) = 1,    𝜗(𝑡1) = 1,   𝜗(𝑡2) = 1,    𝜗(𝑡3) = 1, 

And 𝜏 =  [𝜏1, 𝜏2, … , 𝜏𝑚]3 ∈ 𝑅𝑇, we have  

𝜗(𝜏) =
1

𝜌(𝜏)(𝜌(𝜏) − 1)(𝜌(𝜏) − 2)
 (𝜗(𝜏1),.  .  . 𝜗(𝜏𝑚)). 

Proposition 3.2.1 [22,23]: The density 𝜎(𝜏) is the non negative integer factors defined on trees 

𝑅𝑇 , ∀𝜏 ∈ 𝑅𝑇 satisfy 

1- 𝜎(𝑡1) = 1, 𝜎(𝑡2) = 2 , 𝜎(𝑡3) = 6, 

2- with 𝜏 =  [𝜏1, 𝜏2, … , 𝜏𝑚]3, we have 𝜎(𝜏) =  𝜌(𝜏)(𝜌(𝜏) − 1)(𝜌(𝜏) −

2)(𝜎(𝜏1), … , 𝜎(𝜏𝑚)). 

Proposition 3.2.2 [22,23]: The non-negative integer  휀(𝜏), ∀𝜏 ∈ 𝑅𝑇 satisfy 

1- 𝜖(𝑡1) = 1,   휀(𝑡2) = 1 ,    𝜖(𝑡3) = 1, 

2- For the tree 𝜏 = [𝜏1
𝑣1 , … , 𝜏𝑚

𝑣𝑚]
3

 ∈ 𝑅𝑇,  with distinct 𝜏𝑖 we have 휀(𝜏) = (𝜌(𝜏) −

3)! ∏
1

𝑣𝑖
(

𝜀(𝜏𝑖)

𝜌(𝜏𝑖)!)
𝑣𝑖

 ,𝑚
𝑖=1  where  𝑣𝑖 count similar tree of 𝜏𝑖, 𝑖 = 1, . . . , 𝑚. 

Therefore we can represent B-series (41) as follows: 

𝐵(𝛿, 𝛼) = 𝛿(∅) + ∑
ℎ𝜌(𝜏)

𝜌(𝜏)!
𝛿(𝜏) 휀(𝜏) 𝜎(𝜏) 𝜇(𝜏)(𝛼, 𝛼′, 𝛼′′).𝜏∈𝑅𝑇                                                  (42) 

The previous analysis results in the following theorem. 

Theorem 3.1 [22,23]: The exact solution of (18) is a B-series 

𝛼(𝑥0 + ℎ) = ∑
ℎ𝜌(𝜏)

𝜌(𝜏)!𝜏∈𝑅𝑇 휀(𝜏)𝕌(𝜏)(𝛼0, 𝛼0
′ , 𝛼0

′′),                                                                         (43) 

And the first and second derivatives have the following B-series respectively, 

𝛼′(𝑥0 + ℎ) = ∑
ℎ𝜌(𝜏)−1

(𝜌(𝜏)−1)!𝜏∈𝑅𝑇 휀(𝜏)𝕌(𝜏)(𝛼0, 𝛼0
′ , 𝛼0

′′),                                                                  (44) 

𝛼′′(𝑥0 + ℎ) = ∑
ℎ𝜌(𝜏)−2

(𝜌(𝜏)−2)!𝜏∈𝑅𝑇/[𝑡1] 휀(𝜏)𝕌(𝜏)(𝛼0, 𝛼0
′ , 𝛼0

′′) .                                                          (45) 

Lemma 3.3 [22,23]: We can calculate the function 𝜂𝑖(𝜏) on ∈ 𝑅𝑇/[𝑡1 , 𝑡2 ] recursively as 

1- 𝜂𝑖(𝑡3) = 1, 

2- For the tree  𝜏 = [𝜏1
𝑣1 , 𝑡2

𝑣2 , 𝑡3
𝑣3 , .  .  . , 𝜏𝑚

𝑣𝑚]3  ∈ 𝑅𝑇, with distinct 𝜏𝑖 , 𝑖 = 1, . . . , 𝑚 and 

different 𝑡1,  and 𝑡2, 𝜂𝑗 =
1

2𝑣2
𝑐𝑗

𝑣1+2𝑣2 ∏ (∑ �̂�𝑗𝑘𝜂𝑘(𝜏𝑖))𝑣𝑖𝑠
𝑘=1

𝑚
𝑖=3 . 

Now, we denote the vector 𝜂(𝜏) = (𝜂1(𝜏)(𝜏),.  .  . ,  𝜂𝑠(𝜏)(𝜏))𝑇, ∀𝛼 ∈ 𝑅𝑇\ {𝑡1, 𝑡2 }.  the initial 

weight associated to 𝛼𝑛+1   is indicated by 𝜕(𝜏) and is defined as follows: 

𝜕(𝜏) = ∑ 𝑠𝑖 𝜂𝑖(𝜏)𝑠
𝑖=1 = 𝑠𝑇 𝜂(𝜏), 

𝜕′(𝜏) is indicated to the initial weight associated with 𝛼′
𝑛+1 and is defined as follows: 

𝜕′(𝜏) = ∑ 𝑠𝑖′  𝜂𝑖(𝜏)𝑠
𝑖=1 = 𝑠′𝑇

 𝜂(𝜏),  

and the initial weight associated with 𝛼′′
𝑛+1 indicated by 𝜕′′(𝜏)  and is defined as follows:  



IHJPAS. 37 (2) 2024 

 

395 
 

𝜕′′(𝜏) = ∑ 𝑠𝑖′′  𝜂𝑖(𝜏)𝑠
𝑖=1 = 𝑠′′𝑇

 𝜂(𝜏). 

As a result, we obtain the following fundamental theorem for the numerical solution and RKTDIO 

technique numerical derivatives. 

Theorem 3.2 [22,23]: When we apply the RKTDIO method (39) on the autonomous problem (18) 

yields the numerical solution 𝛼𝑛+1 and numerical derivatives 𝛼′𝑛+1, 𝛼′′𝑛+1 which have the B-

series as follows: 

𝛼𝑛+1 = 𝛼𝑛 + ℎ𝛼𝑛
′ +

1

2
ℎ2𝛼𝑛

′′ + ∑
ℎ𝜌(𝜏)

𝜌(𝜏)!𝜏∈
𝑅𝑇

{𝑡1,𝑡2}

휀(𝜏)𝜎(𝜏)𝜃(𝜏). 𝕌(𝜏)(𝛼0, 𝛼0
′ , 𝛼0

′′),                       (45) 

𝛼𝑛+1
′ = 𝛼𝑛

′ + ℎ𝛼𝑛
′′ + ∑

ℎ𝜌(𝜏)−1)

𝜌(𝜏)!𝜏∈𝑅𝑇/{𝑡1,𝑡2} 휀(𝜏)𝜎(𝜏)𝜃′(𝜏). 𝕌(𝜏)(𝛼0, 𝛼0
′ , 𝛼0

′′),                             (46) 

𝛼𝑛+1
′′ = 𝛼𝑛

′′ + ∑
ℎ𝜌(𝜏)−2)

𝜌(𝜏)!𝜏∈𝑅𝑇/{𝑡1,𝑡2} 휀(𝜏)𝜎(𝜏)𝜃′′(𝜏). 𝕌(𝜏)(𝛼0, 𝛼0
′ , 𝛼0

′′) .                                     (47) 

 

4.  Algebraic order conditions 

 We arrive at this paper's main contribution—the order conditions of the RKTDIO method-through 

Theorems 3.1 and 3.2. In Table 1, the relevant-colored trees of orders up to six are listed together 

with the accompanying function values. 

Theorem 4.1: The RKTDIO method has order 𝑝 (𝑝 ≥ 3) if and only if it satisfies the following 

conditions. 

1- 𝜃(𝜏) =
1

𝜎(𝜏)
, 𝜏 ∈ ⋃ 𝑠𝑟

𝑝
𝑟=4 ,                                                                                               (48) 

2- 𝜃′(𝜏) =
𝜌(𝜏)

𝜎(𝜏)
, 𝜏 ∈ ⋃ 𝑠𝑟

𝑝+1
𝑟=4 ,                                                                                              (49) 

3- 𝜃′′(𝜏) =
𝜌(𝜏)(𝜌(𝜏)−1)

𝜎(𝜏)
, 𝜏 ∈ ⋃ 𝑠𝑟

𝑝+2
𝑟=4 .                                                                                  (50) 

Even though some of the trees in the set RT provide the same order criteria and pertain to various 

elementary differentials, it is still unnecessary. In general, the following corollary, which may be 

derived from the definition of density and order, and from Lemma 3.3, can be used to overcome 

the similarity between the order criteria. 

 

Table 1. lists elementary differentials, relevant-colored trees of up to six orders, and related functions. 

 

Order 𝝉 Tree 𝜶(𝝉) Density 𝒏(𝝉) Elementary 

0 ∅ ∅ 1 1  𝛼 

1 𝑡1  1 1  𝛼′ 

2 𝑡2 
 

1 2  𝛼′′ 

3 𝑡3 

 

1 6  𝜇 

4 𝜏41 

 

1 24 𝑐 𝜇′𝛼′ 



IHJPAS. 37 (2) 2024 

 

396 
 

5 𝜏51 

 

1 60 𝑐2 𝜇′𝛼′′ 

5 𝜏52 

 

1 120 
1

2
𝑐2 𝜇′′(𝛼′, 𝛼′) 

6 𝜏61 

 

3 120 𝑐3 𝜇′′(𝛼′, 𝛼′′) 

6 𝜏62 

 

1 240 
1

2
𝑐3 𝜇′′′(𝛼′, 𝛼′, 𝛼′) 

6 𝜏63 

 

1 720 
1

6
𝑐3 𝜇′𝛼′′ 

 

The order criteria for the RKTDIO technique up to the sixth order can be expressed as follows 

using Theorem 4.1 as a foundation: 

Order condition for 𝜶 

order 3   

∑ 𝑘𝑖 =
1

6
                                                                                                                                     (51) 

Order 4     

∑ 𝑘𝑖 𝑠𝑖 =
1

24
                                                                                                                                (52) 

Order 5   

∑ 𝑘𝑖𝑠𝑖
2 =

1

60
                                                                                                                               (53) 

Order 6 

∑ 𝑘𝑖𝑠𝑖
3 =

1

120
,                      ∑ 𝑘𝑖𝑎𝑖𝑗 =

1

720
                                                                                   (54) 

Order condition for 𝜶′  

Order 2  

∑ 𝑘′𝑖 =
1

2
                                                                                                                                    (55) 

Order 3       

∑ 𝑘𝑖′ 𝑠𝑖 =
1

6
                                                                                                                                 (56) 

Order 4 



IHJPAS. 37 (2) 2024 

 

397 
 

∑ 𝑘𝑖′𝑠𝑖
2 =

1

12
                                                                                                                             (57) 

Order 5 

∑ 𝑘𝑖′𝑠𝑖
3 =

1

20
,                      ∑ 𝑘𝑖′𝑎𝑖𝑗 =

1

120
                                                                                 (58) 

Order 6 

∑ 𝑘𝑖′𝑠𝑖
4 =

1

30
,                     ∑ 𝑘𝑖′𝑎𝑖𝑗𝑠𝑗 =

1

720
,  ∑ 𝑘𝑖′𝑠𝑗𝑎𝑖𝑗 =

1

180
                                                  (59) 

 Order condition for 𝜶′′   

Order 1       

∑ 𝑘𝑖
′′ = 1                                                                                                                                 (60) 

Order 2       

∑ 𝑘𝑖′′ 𝑠𝑖 =
1

2
                                                                                                                              (61) 

Order 3 

∑ 𝑘𝑖′′𝑠𝑖
2 =

1

3
                                                                                                                             (62) 

Order 4 

∑ 𝑘𝑖′′𝑠𝑖
3 =

1

4
,                      ∑ 𝑘𝑖′′𝑎𝑖𝑗 =

1

24
                                                                                  (63) 

Order 5 

∑ 𝑘𝑖′′𝑠𝑖
4 =

1

5
,        ∑ 𝑘𝑖′′𝑎𝑖𝑗𝑠𝑗 =

1

120
,     ∑ 𝑘𝑖′′𝑠𝑗𝑎𝑖𝑗 =

1

30
                                                           (64) 

Order 6 

∑ 𝑘𝑖
′′𝑠𝑖

2𝑎𝑖𝑗 =
1

36
,      ∑ 𝑘𝑖′′𝑠𝑖

5 =
1

6
  , ∑ 𝑘𝑖′′𝑎𝑖𝑗𝑠𝑗

2 =
1

360
,      ∑ 𝑘𝑖′′𝑠𝑖𝑎𝑖𝑗𝑠𝑗 =

1

144
                      (65)  

 

5.  The Construction of RKTDIO Method  

When creating implicit RKTDIO methods, the order conditions listed in Section 3.2 must be met. 

For the 𝑞-order RKTDIO method, the local truncated error is defined as follows: 

‖𝐿𝑔
(𝑞+1)‖2 = (∑ (𝐿𝑖

(𝑞+1) )2 + (∑ (𝐿𝑖′
(𝑞+1) )2 + (∑ (𝐿𝑖′′

(𝑞+1) )2𝑛′′
𝑞+1

𝑖=1

𝑛′
𝑞+1

𝑖=1

𝑛𝑞+1

𝑖=1
                        (66) 

Where 𝐿(𝑞+1) , 𝐿′(𝑞+1), 𝐿′′(𝑞+1) the local truncation error are terms respectively,  𝐿𝑔
(𝑞+1) is the 

global local truncation error.  

5.1   A Three-Stage Fifth-Order RKTDIO Method 

 In this subsection, the derivation of the three-stage RKTDIO technique of order five by using the 

algebraic order conditions up to order five will be considered. The resulting system consists of 16 

nonlinear equations with 16 unknown variables, solving the system simultaneously, and assuming 

𝑎11 = 𝑎22 and 𝑎22 = 𝑎33 

𝑎21 = 𝑎21, 𝑎31 = 𝑎21, 𝑎32 =
3

20
 𝑅𝑜𝑜𝑡 𝑂𝑓 (5 − 𝑧2 − 3), 𝑎33 = −

5

9
 𝑎21 −

1

24
 𝑅𝑜𝑜𝑡 𝑂𝑓(5 − 𝑧2 −

3) +
1

24
 , 𝑏1 =

2

9
 , 𝑏2 =

1

36
 
200 𝑎21𝑅𝑜𝑜𝑡 𝑂𝑓(5−𝑧2−3)+15𝑅𝑜𝑜𝑡 𝑂𝑓(5−𝑧2−3)+200 𝑎21+9 

40 𝑎21+3 𝑅𝑜𝑜𝑡 𝑂𝑓( 5−𝑧2−3)
, 𝑏3 =

−
1

36
 
200 𝑎21𝑅𝑜𝑜𝑡 𝑂𝑓(5−𝑧2−3)−15𝑅𝑜𝑜𝑡 𝑂𝑓(5−𝑧2−3)+200 𝑎21+9 

40 𝑎21+3 𝑅𝑜𝑜𝑡 𝑂𝑓( 5−𝑧2−3)
, 𝑐1 =

1

2
, 𝑐2 =

1

2
−

1

2
 𝑅𝑜𝑜𝑡 𝑂𝑓( 5 − 𝑧2 −

3), 𝑐3 =
1

2
 𝑅𝑜𝑜𝑡 𝑂𝑓( 5 − 𝑧2 − 3) +

1

2
, 𝑑1 =

1

18
, 𝑑2 =  

5

72
 𝑅𝑜𝑜𝑡 𝑂𝑓( 5 − 𝑧2 − 3) +

1

18
, 𝑑3 =

− 
5

72
 𝑅𝑜𝑜𝑡 𝑂𝑓( 5 − 𝑧2 − 3) +

1

18
, 𝑔1 =

4

9
, 𝑔2 =

5

18
, 𝑔3 =

5

18
. 

Next, we minimize the truncation error term by using minimize command in Maple. Thus, for the 

optimized value of coefficients in fractional we chose  𝑎21 = −
1

125
 with this value ‖𝜋𝑔

(5)‖
2

=



IHJPAS. 37 (2) 2024 

 

398 
 

0.001366343866. Finally, all the parameters of three-stage fifth-order RKTDIO approach that 

will be denoted as RKTDIO5 can be written as follows (see Table 2): 

 

Table 2. The RKTDIO5 Method 

 

1

2
 

83

1800
−

√15

120
 

  

1

2
−

√15

10
 −

1

125
 

83

1800
−

√15

120
 

 

1

2
+

√15

10
 −

1

125
 

3√15

100
 

83

1800
−

√15

120
 

 1

18
 

1

18
+

√15

72
 

1

18
−

√15

72
 

 2

9
 

5

36
+

√15

36
 

5

36
−

√15

36
 

 4

9
 

5

18
 

5

18
 

 

5.2   A Four-Stage RDTDIO Method of  Order Six 

        For the four-stage RKTDIO technique of order six, the algebraic conditions up to order six 

will be solved. The resulting system consists of 25 nonlinear equations with 23 unknown variables, 

and using simplifying assumption 𝑏𝑖
′ = 𝑏𝑖

′′(1 − 𝑐𝑖), 𝑖 = 1, … , 𝑠,  and supposing  𝑏2
′′ = 0, and 𝑏2

′′′ =

0  

𝑎21 =
3

20
 𝑅𝑜𝑜𝑡 𝑂𝑓 (10 − 𝑧2 − 10 − 𝑧 + 1) −

3

80
, 𝑎31 = −

1

40
, 𝑎32 =

11

80
−

3

20
 𝑅𝑜𝑜𝑡 𝑂𝑓(10 −

𝑧2 − 10 − 𝑧 + 1), 𝑎41 =
1

40
, 𝑎42 = −

1

40
, 𝑎43 =

3

20
 𝑅𝑜𝑜𝑡 𝑂𝑓(10 − 𝑧2 − 10 − 𝑧 + 1) −

3

80
,

𝑎44 =
1

48
, 𝑏1 =

2

9
, 𝑏2 = 0, 𝑏3 =

5

18
  𝑅𝑜𝑜𝑡 𝑂𝑓(10 − 𝑧2 − 10 − 𝑧 + 1), 𝑏4 = −

5

18
  𝑅𝑜𝑜𝑡 𝑂𝑓(10 −

𝑧2 − 10 − 𝑧 + 1) +
5

36
, 𝑐1 =

1

2
, 𝑐2 = 𝑅𝑜𝑜𝑡 𝑂𝑓(10 − 𝑧2 − 10 − 𝑧 + 1), 𝑐3 = 1 − 𝑅𝑜𝑜𝑡 𝑂𝑓(10 −

𝑧2 − 10 − 𝑧 + 1), 𝑐4 = 𝑅𝑜𝑜𝑡 𝑂𝑓(10 − 𝑧2 − 10 − 𝑧 + 1), 𝑑1 =
1

18
, 𝑑2 = 𝑑2, 𝑑3 =

5

36
𝑅𝑜𝑜𝑡 𝑂𝑓(10 − 𝑧2 − 10 − 𝑧 + 1) −

1

72
, 𝑑4 = −𝑑2 −

5

36
𝑅𝑜𝑜𝑡 𝑂𝑓(10 − 𝑧2 − 10 − 𝑧 + 1) −

1

8
, 𝑔1 =

4

9
, 𝑔2 = 0, 𝑔3 =

5

18
, 𝑔4 =

5

18
  

Lastly, all the parameters of four-stage sixth-order RKTDIO method indicated by RKTDIO6 can 

be written as follows: 

 

Table 3. The RKTDIO6 Method. 

1

2
 

1

48
    

1

2
−

√15

10
 

3

80
−

3√15

200
 

1

48
   

1

2
+

√15

10
 −

1

40
 

1

16
+

3√15

200
 

1

48
  

1

2
−

√15

10
 

1

40
 −

1

40
 

3

80
−

3√15

200
 

1

18
+

√15

72
 



IHJPAS. 37 (2) 2024 

 

399 
 

 
1

18
 0 

1

18
−

√15

72
 

1

18
−

√15

72
 

 
2

9
 0 

5

36
−

√15

36
 

5

36
+

√15

36
 

 
4

9
 0 

5

18
 

5

18
 

 

 

6.  Derivation Embedded ERKTDIO6(5) Method 

        The RKTDIO method with 𝔰-stages for solving equation (1) is presented in its general form. 

Subsequently, the development of the embedded pair RK approach is discussed, which is an area 

of active research aimed at improving existing codes. To estimate the minimum error, implicit 

RKTDIO techniques are used, employing pairs of 𝒫(𝒬) orders in the values of step size codes. 

These methods are based on the 𝒫-order method (𝐶, 𝐴, 𝑏, 𝑏′, 𝑏′′) and the 𝒬-order method 

(𝐶, 𝐴, �̈�, �̈�′, �̈�′′), and can be represented using the Butcher Tabular notation. The embedded pair 

can be initialized in the following manner: 

C A 
  

𝑏 𝑇  
𝑏′ 𝑇  
𝑏′′ 𝑇  

  

�̈� 𝑇  

�̈�′ 𝑇  

�̈�′′ 𝑇  
 
 

 

The main objective of constructing the embedded pair of implicit RKTDIO techniques is to 

obtain a single error estimate that can be used in step-size approaches. This is achieved by 

improving the existing pairs and local error estimates and then restricting the step size 𝓀. 

𝒽𝑛+1 = 0.9𝓀𝑛(
𝑇𝑜𝑙

𝐿𝑇𝐸
 )

1

𝒬+1                                                                                                            (67) 

A safety factor of 0.9 is used to determine the local error estimate at each step, and 𝑇𝑜𝑙 represents 

the maximum allowable local error that ensures the necessary accuracy. If the local truncation 

error (LTE) is less than or equal to 𝑇𝑜𝑙, the step is accepted and the higher order method (or local 

extrapolation) is used, where a more accurate approximation is applied to drive the integration 

and update 𝓀 using Equation (67). On the other hand, if LTE is greater than 𝑇𝑜𝑙, the step is 

rejected and the step size 𝓀 is reduced by half. The RKTDIO method is an embedded Runge-

Kutta method developed for solving third-order ODEs. To ensure high accuracy for the higher-

order method and the most accurate error estimates for the lower-order methods, fractions were 

used to develop orders 6 and 5, respectively. The step size 𝓀 plays a crucial role in obtaining 

accurate results and can be doubled to achieve this goal. The embedded RKTDIO6(5) method is 

derived in Table 4 for this study. 

In RKTDIO6(5), the 𝐴 and 𝐶 values is computed from the 6th-order solution then derived the 

four-stage 5rd-order embedded equation. Solving of the eqs. (51-53), (55-58), and (60-64) 



IHJPAS. 37 (2) 2024 

 

400 
 

simultaneously then the solution for �̈� and �̈�′ while �̈�′′ have the same values as the 6th -order. 

The solutions are obtained as 

�̈�1
′ =

2

9
,  �̈�2

′ = 0,  �̈�3
′ =

5

36
−

√15

36
,  �̈�4

′ =
5

36
+

√15

36
, �̈�1 =

1

18
, �̈�2 =

1

18
+

√15

72
, �̈�3 =

1

18
−

√15

72
, �̈�4 =

�̈�4, �̈�1
′′ =

4

9
, �̈�2

′′ = 0, �̈�3
′′ =

5

18
, �̈�4

′′ =
5

18
 .  

The following simplifying assumption is used in order to reduce the number of equations to be 

solved: 

𝑏𝑖
′ = 𝑏𝑖

′′(1 − 𝑐𝑖), 𝑖 = 1, … , 𝔰.                                                                                                   (68) 

Initially, we have 16 nonlinear equations with 12 unknown variables that we need to find a 

solution for. However, as the number of equations is greater than the number of unknowns, there 

is no solution. To tackle this issue, we make the simplifying assumption (68) which reduces the 

number of equations to 12 with 11 unknowns, enabling us to solve the system. We choose �̈� =
1

10
 as the free parameter, and as a result, we can express the coefficients of the 4-stage embedded 

ERKTDIO6(5) technique. 

 

Table 4: Table of ERKTDIO6(5) method. 

 
7. Numerical Experiments 

7.1 Constant Methods   

     In order to evaluate the performance of the new RKTDIO methods to the established RK 

methods in the scientific literature, a series of test problems are addressed in this part. The 

following methods have been selected for comparison: 

 RKTDIO5: the implicit RKTDIO method of order five with three-stage derived in this 

paper. 

 RKTDIO6: the implicit RKTDIO method of order six with four-stage derived in this 

paper. 



IHJPAS. 37 (2) 2024 

 

401 
 

 DITRKM5: three-stage fifth-order implicit RK derived by [25]. 

 Radau I: three-stage fifth-order implicit RK method presented in [28]. 

 Radau IA: three-stage fifth-order implicit RK method presented in [24]. 

 Lobatto III: The sixth-order four-stage implicit Runge–Kutta method as given by 

[30]. 

 Lobatto IIIB: The sixth-order four-stage implicit Runge–Kutta method as given by 

[29]. 

Problem (1): Consider a nonhomogeneous linear ODE given in [22] 

α′′′(x) = α(x) + cos(x),  

With  

 𝛼(0) = 0 , 𝛼′(0) = 0, 𝛼′′(0) = 1 

Where x ∈ [0,1],  

  and analytic solution   𝛼(𝑥) =
(𝑒𝑥−cos(𝑥)−sin(𝑥))

2
. 

 

Problem (2): Consider the nonhomogeneous nonlinear ODE 

α′′′(x) = (α(x))2 + cos2(x) − cos(x) − 1 , 

With 

𝛼(0) = 0, 𝛼′(0) = 1, 𝛼′′(0) = 1  𝑤ℎ𝑒𝑟𝑒 0 ≤ 𝑥 ≤ 2 , 

and analytic solution 𝛼(𝑥) = sin (𝑥). 

 

Problem (3): The nonhomogeneous nonlinear ODEs is considered as 

α′′′(x) = 8(
𝛼2(𝑥)

𝑒2𝑥
) 

With 

𝛼(0) = 1, 𝛼′(0) = 2 , 𝛼′′(0) = 4    𝑤ℎ𝑒𝑟𝑒 0 ≤ 𝑥 ≤ 1 , 

and analytic solution  𝛼(𝑥) = 𝑒2𝑥 . 

 
Figure 2. Accuracy curve for RKTDIO5, DITRKM5, Radau I, Radau IA with ℎ =  0.1, 0.05, 0.025, 0.00125, 0.00625 

for Problem1. 



IHJPAS. 37 (2) 2024 

 

402 
 

 
Figure3. Accuracy curve for RKTDIO5, DITRKM5, Radau I, Radau IA with ℎ =  0.1, 0.05, 0.025, 0.00125, 0.00625 

for Problem2. 

 

 
Figure4. Accuracy curve for RKTDIO5, DITRKM5, Radau I, Radau IA with ℎ =  0.1, 0.05, 0.025, 0.00125, 0.00625 

for Problem3. 

 
Figure 5. Accuracy curve for RKTDIO6, Lobattoo III, Lobattoo IIIB with ℎ =  0.1, 0.05, 0.025, 0.00125, 0.00625 

for Problem1. 



IHJPAS. 37 (2) 2024 

 

403 
 

 
Figure 6. Accuracy curve for RKTDIO6, Lobattoo III, Lobattoo IIIB with ℎ =  0.1, 0.05, 0.025, 0.00125, 0.00625 

for Problem2. 

 
Figure 7. Accuracy curve for RKTDIO6, Lobattoo III, Lobattoo IIIB with ℎ =  0.1, 0.05, 0.025, 0.00125, 0.00625 

for Problem3. 

7.2 Variable Method 

         This subsection will apply the new embedded method to the same third-order differential 

equation problems in the previous subsection. The following implicit diagonally RK method is 

selected for the numerical comparisons.  The approximation results are illustrated in the tables 

below for solving problems. The following abbreviations will be used in the tables: 

 𝑻𝒐𝒍: Tolerance. 

 Method: method employed step sizes between two points or positions. 

 F. N: number of the function call. 

 STEP: The number of successful steps. 

 FSTEP: The number of failed steps. 

 Time: execution time. 

 ERKTDIO6(5): The novel embedded 6(5) derived in this study. 

 EDITRK5(4): The embedded diagonally implicit 5(4) RK derived in [7]. 

 



IHJPAS. 37 (2) 2024 

 

404 
 

Table 4.  Comparisons of number of function call and Time of ERKTDIO6(5) and EDITRK 5(4) with  h =

10−6, 10−8, 10−10 for the problem 1. 

FSTEP Step Time No. of Function Call Method 𝑻𝑶𝑳(𝓴) 

2 5 0.025 270 ERKTDIO6(5) 
10−6 

0 13 0.033  421 EDITRK 5(4) 

2 16 0.038 1210 ERKTDIO6(5) 
10−8 

0 40 0.047  1962 EDITRK 5(4) 

2 72 0.060 6353 ERKTDIO6(5) 
10−10 

1  128 0.074 9105 EDITRK 5(4) 

 

Table 5. Comparisons of number of function call and Time of EDITRKM 4(3) and EDITRKM 5(4) with h =

10−2, 10−4, 10−6 for the problem 2. 

FSTEP             Step Time No. of Function Call            Method 𝑻𝑶𝑳(𝓴) 

0 3 0.032 183 ERKTDIO6(5) 
10−6 

0 7 0.062 166 EDITRK 5(4) 

0 8 0.045 662 ERKTDIO6(5) 
10−8 

0 22 0.083 759 EDITRK 5(4) 

0 29 0.080 2533 ERKTDIO6(5) 
10−10 

1 69 0.119 3495 EDITRK 5(4) 

 

Table 6. Comparisons of number of function call and Time of EDITRKM 4(3) and EDITRKM 5(4) with h =

10−2, 10−4, 10−6 for the problem 3. 

FSTEP             Step Time No. of Function Call            Method 𝑻𝑶𝑳(𝓴) 

1 97 0.012 466 ERKTDIO6(5) 
10−6 

0 29 0.035 882 EDITRK 5(4) 

1 315 0.031 2511 ERKTDIO6(5) 
10−8 

1 152 0.059 4101 EDITRK 5(4) 

1 1260 0.075 13325 ERKTDIO6(5) 
10−10 

2 725 0.095 19044 EDITRK 5(4) 

 

 

 



IHJPAS. 37 (2) 2024 

 

405 
 

 

Problem 1 

 

Problem 2 

 

Problem 3 

Figure 8: Accuracy curve for ERKTDIO6(5) and EDITRK5(4) with ℎ = 10−6, 10−8, 10−10 . 



IHJPAS. 37 (2) 2024 

 

406 
 

8. Discussions  

      The objective of our research was to directly solve third-order ordinary differential equations 

(ODEs) using the RKTDIO approach, which is based on the algebraic theory of order conditions. 

Previous research has primarily concentrated on defining algebraic order conditions and B-series 

theory to solve first- and second-order ordinary differential equations (ODEs). However, we were 

inspired to expand upon these ideas and create the RKTDIO formula specifically designed for 

third-order ODEs. As a result, the RKTDIO5 and RKTDIO6 methods were developed, offering 

improved computing efficiency and accuracy for solving certain third-order ODEs as compared to 

other current approaches. 

      In addition, we shared the results of our study on the embedded diagonal implicit type Runge-

Kutta technique (ERKTDIO). We assessed the performance of our method (ERKTDIO6(5)) in 

comparison to other approaches by examining the decimal logarithm of the highest time curve and 

the logarithm of the function call estimates obtained from Tables 4-6. We utilised three separate 

test problems and calculated the logarithm of the time curve utilising various tolerance values. 

Figure 10 was generated using the numerical data obtained from Tables 4-6. It presents the number 

of successful and failed steps that occurred during the calculations. Our work has enhanced and 

refined the method by transitioning from an explicit to an implicit approach and from a direct to a 

diagonal scheme. 

 

9. Conclusions 

 Ultimately, our research has made a significant contribution to the progress of numerical methods 

for solving third-order ordinary differential equations. The introduction of the RKTDIO5 and 

RKTDIO6 techniques provides improved computational efficiency and precision for particular 

third-order ordinary differential equations (ODEs). Furthermore, our examination of the 

ERKTDIO technique showcases its efficacy in comparison to alternative methods, especially in 

terms of the number of successful steps in computations. Our methods provide the potential for 

enhanced solutions of third-order ODEs by adopting implicit and diagonal schemes. In summary, 

these discoveries expand the current understanding of numerical analysis and computational 

mathematics, offering useful insights for future investigations in this field. 

 

Acknowledgment 

The authors have 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. Boutayeb, A.; Chetouani A. A mini-review of numerical methods for high-order problems. International 

Journal of Computer Mathematics. 2007, 84(4), 563-579, https://doi.org/10.1080/00207160701242250. 

https://doi.org/10.1080/00207160701242250


IHJPAS. 37 (2) 2024 

 

407 
 

2. Wu, Xiong-Jian; Yigong Wang; W. G. Price. Multiple resonances, responses, and parametric 

instabilities in offshore structures. Journal of ship research. 1988, 32(04), 285-296, 

https://doi.org/10.5957/jsr.1988.32.4.285 

3. Cortell, R. Application of the fourth-order Runge-Kutta method for the solution of high-order general 

initial value problems. Computers & structures. 1993, 49(5), 897-900, https://doi.org/10.1016/0045-

7949(93)90036-D 

4. Malek, A.; R. Shekari Beidokhti. Numerical solution for high order differential equations using a hybrid 

neural network-optimization method.  Applied Mathematics and Computation. 2006, 183(1), 260-271, 

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

5. Awoyemi, D. O.  Algorithmic collocation approach for direct solution of fourth-order initial-value 

problems of ordinary differential equations.  International Journal of Computer Mathematics. 

2005, 82(3), 321-329, https://doi.org/10.1080/00207160412331296634. 

6. Adel, F.; Senu, N.; Ismail, F.; Majid, Z. A. New efficient phase-fitted and amplification-fitted Runge-

Kutta method for oscillatory problems. International J Pure Appl Math, 2016, 107(1), 69-86, 

https://doi.org/10.12732/ijpam.v107i1.7 

7. Hussain, K.; Ismail, F., ; Senu, N. A new optimized Runge-Kutta-Nyström method to solve oscillation 

problems. World Applied Sciences Journal, 2015, 33(10), 1614-1622, 

https://doi.org/10.5829/idosi.wasj.2015.33.10.300. 

8. Hussain, K. A.; Abdulnaby, Z. E.. A new two derivative FSAL Runge-Kutta method of order five in 

four stages. Order, 2020, 17(1). https://doi.org/10.21123/bsj.2020.17.1.0166 

9. Lee, K. C.; Senu, N., Ahmadian, A., Ibrahim, S. N. I., & Baleanu, D. Numerical study of third-order 

ordinary differential equations using a new class of two derivative Runge-Kutta type 

methods. Alexandria Engineering Journal, 2020, 59(4), 2449-2467, 

https://doi.org/10.1016/j.aej.2020.03.008. 

10. Ahmad, N. A.; Senu, N.; Ismail, F.  Phase-fitted and amplification-fitted higher order two-derivative 

Runge-Kutta method for the numerical solution of orbital and related periodical IVPs. Mathematical 

Problems in Engineering, 2017. 10.1109/ICREM.2015.7357020. 

11. Senu, N.; Ahmad, N. A., ; zarina, B. I. B. I. Numerical Study on Phase-Fitted and Amplification-Fitted 

Diagonally Implicit Two Derivative Runge-Kutta Method for Periodic IVPs. Sains Malaysiana, 2021, 

50(6), 1799-1814, 10.17576/jsm-2021-5006-25. 

12. Fawzi, F. A. An Embedded 5 (4) Pair of Optimized Runge-Kutta Method for the Numerical Solution of 

Periodic Initial Value Problems. Iraqi Journal of Science, 2022,  3889-3901, 

https://doi.org/10.24996/ijs.2022.63.9.21. 

13. SENU, N.; AHMAD, N.;  ISMAIL, F. Embedded 5 (4) Pair Trigonometrically-Fitted Two Derivative 

Runge-Kutta Method with FSAL Property for Numerical Solution of Oscillatory Problems, WSEAS 

TRANSACTIONS on COMPUTERS, 2017,16, 155-162. 

14. Jikantoro, Y. D.; Ismail, F., Senu, N., ; Ibrahim, Z. B.,  Hybrid methods for direct integration of special 

third order ordinary differential equations. Applied Mathematics and Computation, 2018,  320, 452-463, 

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

15. Shokri, A.; Mehdizadeh Khalsaraei, M. A new efficient high order four-step multiderivative method for 

the numerical solution of second-order IVPs with oscillating solutions. Mathematics Interdisciplinary 

Research, 2020, 5(2), 157-172, 10.22052/mir.2020.211603.1185. 

16. Mechee, M.; Senu, N.; Ismail, F.; Nikouravan, B.; Siri, Z. A three-stage fifth-order Runge-Kutta method 

for directly solving special third-order differential equation with application to thin film flow problem. 

Mathematical Problems in Engineering. 2013, 2013,  https://doi.org/10.1155/2013/795397. 

17. Fawzi, F. A.; Jumaa M, H. The Implementations Special Third-Order Ordinary Differential Equations 

(ODE) for 5th-order 3rd-stage Diagonally Implicit Type Runge-Kutta Method (DITRKM). Ibn AL-

https://doi.org/10.5957/jsr.1988.32.4.285
https://doi.org/10.1016/0045-7949(93)90036-D
https://doi.org/10.1016/0045-7949(93)90036-D
https://doi.org/10.1016/j.amc.2006.05.068
https://doi.org/10.1080/00207160412331296634
http://dx.doi.org/10.12732/ijpam.v107i1.7
https://doi.org/10.21123/bsj.2020.17.1.0166
https://doi.org/10.1016/j.aej.2020.03.008
http://dx.doi.org/10.1109/ICREM.2015.7357020
http://dx.doi.org/10.17576/jsm-2021-5006-25
https://doi.org/10.24996/ijs.2022.63.9.21
https://doi.org/10.1016/j.amc.2017.10.003
http://dx.doi.org/10.22052/mir.2020.211603.1185
https://doi.org/10.1155/2013/795397


IHJPAS. 37 (2) 2024 

 

408 
 

Haitham Journal For Pure and Applied Sciences. 2022, 35(1), 92-101, 

https://doi.org/10.30526/35.1.2803.  

18. ISMAIL F.; Al-Khassawneh, R. A. Solving Delay Differential Equations Using embedded singly 

diagonally Implicit runge-kutta methods. Acta Mathematica Vietnamica. 2008, 33(2), 95-105. 

19. Jikantoro, Y. D., Ismail, F., Senu, N., & Ibrahim, Z. B,. A new integrator for special third order 

differential equations with application to thin film flow problem. Indian Journal of Pure and Applied 

Mathematics, 2018, 49, 151-167, 10.1007/s13226-018-0259-6. 

20. Lee, K. C.; Senu, ; Ahmadian, A. ; Ibrahim, S. N. I., On two-derivative Runge–Kutta type methods for 

solving u‴= f (x, u (x)) with application to thin film flow problem. Symmetry, 2020, 12(6), 

924,  https://doi.org/10.3390/sym12060924. 

21. Hussain, K.; Ismail, F.; Senu, N. Solving directly special fourth-order ordinary differential equations 

using Runge–Kutta type method. Journal of Computational and Applied Mathematics. 2016, 306, 179-

199, https://doi.org/10.1016/j.cam.2016.04.002. 

22. Ghawadri, N.;  Senu, N.; Adel Fawzi, F.; Ismail, F.; Ibrahim, Z. B. Explicit Integrator of Runge-Kutta 

Type for Direct Solution of 𝑢(4) =  𝑓 (𝑥, 𝑢, 𝑢′, 𝑢′′). Symmetry. 2019, 11(2), 246, 

https://doi.org/10.3390/sym11020246.  

23. Dormand, J. R. Numerical methods for differential equations: a computational approach. CRC press, 

1996, https://doi.org/10.1201/9781351075107 . 

24. Fawzi, F. A.; Jaleel, N. W. Exponentially Fitted Diagonally Implicit EDITRK Method for Solving 

ODEs. Ibn AL-Haitham Journal For Pure and Applied Sciences. 2023, 36(1), 389-399, 

https://doi.org/10.30526/36.1.2883 . 

25. Saleh, O. M.; Fawzi, F. A.; Hussain, K. A.; Ghawadri, N. G. Derivation of Embedded Diagonally 

Implicit Methods for Directly Solving Fourth-order ODEs. Ibn AL-Haitham Journal For Pure and 

Applied Sciences, 2024, 37(1), 375-385, https://doi.org/10.30526/37.1.3288. 

26. Senu, N.; Suleiman, M.; Ismail, F. An embedded explicit Runge–Kutta–Nyström method for solving 

oscillatory problems. Physica Scripta. 2009, 80(1), 015005.  https://doi.org/10.1088/0031-

8949/80/01/015005 

27. Hairer E.; Hochbruck, M.; Iserles, A.; Lubich, C. Geometric numerical integration; Oberwolfach 

Reports. 2006, 3(1), 805-882. https://doi.org/10.1017/S0962492902000144. 

28. Hairer E, Norsett SP, Wanner G. Solving Ordinary Differential Equations I: Nonstiff Problems, II: Stiff 

Problems. Springer, 1987. https://doi.org/10.1007/978-3-540-78862-1. 

29. Lambert, J. D. Numerical methods for ordinary differential systems. New York: Wiley, 1991. 

https://doi.org/10.1002/9781119121534 

https://doi.org/10.30526/35.1.2803
http://dx.doi.org/10.1007/s13226-018-0259-6
https://doi.org/10.3390/sym12060924
https://doi.org/10.1016/j.cam.2016.04.002
https://doi.org/10.3390/sym11020246
https://doi.org/10.1201/9781351075107
https://doi.org/10.30526/36.1.2883
https://doi.org/10.30526/37.1.3288
http://dx.doi.org/10.1088/0031-8949/80/01/015005
http://dx.doi.org/10.1088/0031-8949/80/01/015005
http://dx.doi.org/10.1017/S0962492902000144
http://dx.doi.org/10.1007/978-3-540-78862-1