375 

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

 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 

 

 

 

 

   

1Omar Mohammed Saleh    2 Firas Adel Fawzi    3Kasim Abbas Hussain*  

4Nizam G. Ghawadri  

 
1,2Department of Mathematics, Faculty of Computer Science and Mathematics, University of Tikrit, 

Tikrit, Iraq. 
3Department of Mathematics, College of Science, Mustansiriyah University, Baghdad, Iraq. 

                                            4Doctor of Mathematics, Ministry of Education, Jenin, Palestine. 

 

 

*Corresponding Author. kasimabbas@uomustansiriyah.edu.iq 

 

 

 

 

Abstract   

         EDIRKTO is an Implicit Type Runge-Kutta  Method of Diagonally Embedded pairs, is a 

novel approach presented in the paper that may be used to solve 4th-order ordinary differential 

equations of the form 𝑞(4)(𝑡)  =  𝑓(𝑡, 𝑞(𝑡), 𝑞′(𝑡)). There are two pairs of EDIRKTO, with three 

stages each: EDIRKTO4(3) and EDIRKTO5(4). The derivation techniques of the method 

indicate that the higher-order pair is more accurate, while the lower-order pair provides superior 

error estimates. Next, using these pairs as a basis, we developed variable step codes and applied 

them to a series of 4𝑡ℎ-order ODE problems. The numerical outcomes demonstrated how much 

more effective their approach is in reducing the quantity of function evaluations needed to 

resolve fourth-order ODE issues. 

 Keywords: Fourth-order ODEs; Runge-Kutta methods; diagonally implicit technique; 

embedded method. 

 

1. Introduction 

Differential equations (DEs) are important tools in mathematical modeling, particularly in 

explaining physical phenomena like heat and fluid movement, object motion, and nuclear 

reactions. ODEs of different orders are used in various applied research fields such as 

mechanics, electrical and control engineering, fluid dynamics, ship dynamics, neural networks, 

and beam theory [1-3]. Numerical and approximated methods have been developed to solve 

specific types of DEs. 

Received 17 February 2023, Received 15 April 2023, Accepted 10 May 2023, Published 20 January 2024 

 

 

Derivation of Embedded Diagonally Implicit Methods for Directly Solving 

Fourth-order ODEs 

 

doi.org/ 10.30526/37.1.3288 

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
mailto:kasimabbas@uomustansiriyah.edu.iq
https://orcid.org/0009-0003-6837-6231
mailto:Oommaarr7ag7@gmail.com
https://orcid.org/0000-0003-0939-7940
mailto:firasadil01@tu.edu.iq
https://orcid.org/0000-0001-9389-3861
mailto:kasimabbas@uomustansiriyah.edu.iq
https://orcid.org/0000-0002-4335-1789
mailto:nizamghawadri@gmail.com


 IHJPAS. 37 (1)  2024 

376 
 

In this study, we discuss a class of quasi-linear, fourth-order (ODEs) and their numerical 

integration in the following form: 

           𝑞(4)(𝑡) = 𝑓(𝑡, 𝑞(𝑡), 𝑞′(𝑡)),    𝑡 ≥ 𝑡0,                                                                            (1) 

with initial conditions 

𝜎𝑖(𝜏)  =  𝜁𝑖, 𝑖 = 0, 1, … , 3. 

Where 𝑓 ∶  ℛ × ℛ𝑁  ⟶ ℛ𝑁 , 𝜎(𝜏)  =  [𝜎1(𝜏), 𝜎2(𝜏), … , 𝜎𝑁(𝜏)] , 𝑡(𝜏, 𝜎) =  [𝑡1(𝜏, 𝜎), 𝑡2(𝜏, 𝜎),

… , 𝑡𝑁(𝜏, 𝜎)], and  𝜁𝑖 =  [𝜁1
i , 𝜁2

i , … , 𝜁N
i  ], where i = 0,1,2, … ,4. 

     Previously, researchers transformed Equation (1) into a system of first-order ODEs with four 

additional dimensions to solve it. However, it would be more efficient if numerical methods could 

be used to solve the problem accurately and quickly. In [4-12] contain such works. Multistep 

strategies for solving ODEs require initial values. However, according to several researchers (see 

[1, 13, 14]) this technique has a drawback as it consumes a lot of computing time and human effort. 

Consequently, several researchers have turned their attention to direct integration methods for 

solving higher-order ODEs, as these methods have demonstrated features of accuracy and speed 

[16-30]. Nevertheless, all of the methods mentioned above are above are multistep. 

     The primary objective of this paper is to introduce a new technique named EDIRKTO, which 

is a one-step implicit Runge-Kutta method created to directly solve general 4𝑡ℎ-order differential 

equations. The method involves deriving diagonally embedded implicit Runge-Kutta methods for 

the direct integration of specific 4𝑡ℎ-order differential equations. To address the IVPs problem in 

(1), the special version of the EDIRKTO method with 𝑚 stages is as follows:  

𝑞𝑛+1 = 𝑞𝑛 + ℎ 𝑞𝑛
′ +

ℎ2

2
 𝑞𝑛

′′ +
ℎ3

6
 𝑞𝑛

′′′ + ℎ4 ∑ 𝑏𝑖𝑓(𝑡𝑛 + 𝑐𝑗ℎ𝑠
𝑖=1 , 𝑄𝑖, 𝑄𝑖

′),                                  (2) 

𝑞𝑛+1
′ =  𝑞𝑛

′ + ℎ 𝑞𝑛
′′ +

ℎ2

2
 𝑞𝑛

′′′ + ℎ3 ∑ 𝑏𝑖
′ 𝑓(𝑡𝑛 + 𝑐𝑗ℎ, 𝑄𝑖, 𝑄𝑖

′),𝑠
𝑖=1                                               (3) 

𝑞𝑛+1
′′ =  𝑞𝑛

′′ + ℎ 𝑞𝑛
′′′ + ℎ2 ∑ 𝑏′𝑖

′ 𝑓(𝑡𝑛 + 𝑐𝑗ℎ, 𝑄𝑖, 𝑄𝑖
′),𝑠

𝑖=1                                                            (4) 

𝑞𝑛+1
′′′ =  𝑞𝑛

′′′ + ℎ ∑ 𝑏′′𝑖
′ 𝑓(𝑡𝑛 + 𝑐𝑗ℎ, 𝑄𝑖, 𝑄𝑖

′),𝑠
𝑖=1                                                                         (5) 

𝑄𝑗 = 𝑞𝑛 + ℎ 𝑐𝑖𝑞𝑛
′ +

ℎ2

2
 𝑐𝑖

2𝑞𝑛
′′ +

ℎ3

6
 𝑐𝑖

3𝑞𝑛
′′′ + ℎ4 ∑ 𝑎𝑖𝑗𝑓(𝑡𝑛 + 𝑐𝑗ℎ, 𝑄𝑗, 𝑄𝑗

′)𝑠
𝑗=1 ,                         (6) 

𝑄𝑗
′ = 𝑞𝑛

′ + ℎ 𝑐𝑖𝑞𝑛
′′ +

ℎ2

2
 𝑐𝑖

2𝑞𝑛
′′′ + ℎ3 ∑ �̅�𝑖𝑗𝑓(𝑡𝑛 + 𝑐𝑗ℎ, 𝑄𝑗 , 𝑄𝑗

′)𝑠
𝑗=1 .                                          (7) 

The coefficients  𝑏𝑖, 𝑏𝑖
′, 𝑏𝑖

′′, 𝑏𝑖
′′′, 𝑎𝑖,𝑗 , �̅�𝑖𝑗  and  𝑐𝑖 of diagonal implicit RK type (EDIRKTO) 

methods are real numbers. The method is diagonally implicit when  𝑎𝑖,𝑗  ≠ 0  for  𝑖 ≤  𝑗. Using 

the Butcher tableau, the EDIRKT approach is illustrated. See Table 1. 

 

 

 

 

 



IHJPAS. 37 (1) 2024 

377 
 

Table 1. The Butcher tableau EDIRKTO method 

 

3. Order Conditions of the EDIRKTO Method 

According to Ghawadri [15], the algebraic order conditions for the EDIRKTO formula up to order 

seven are as follows: 

order 1:  ∑ 𝑏𝑖
′′′ = 1, 

order 2:  ∑ 𝑏𝑖
′′′𝑐𝑖 =

1

2
,       ∑ 𝑏𝑖

′′ =
1

2
, 

order 3:  ∑ 𝑏𝑖
′′′𝑐𝑖

2 =
1

3
,     ∑ 𝑏𝑖

′′𝑐𝑖 =
1

6
,      ∑ 𝑏𝑖

′ =
1

6
, 

order 4:  ∑ 𝑏𝑖
′′′𝑐𝑖

3 =
1

4
,     ∑ 𝑏𝑖

′′′𝑎𝑖𝑗̅̅ ̅̅ =
1

24
,   ∑ 𝑏𝑖

′′𝑐𝑖
2 =

1

12
,     ∑ 𝑏𝑖

′𝑐𝑖 =
1

24
,    ∑ 𝑏𝑖 =

1

24
, 

order 5:  ∑ 𝑏𝑖
′′′𝑐𝑖

4 =
1

5
,   ∑ 𝑏𝑖

′′′𝑎𝑖𝑗 =
1

120
,    ∑ 𝑏𝑖

′′′�̅�𝑖𝑗𝑐𝑗 =
1

120
,   ∑ 𝑏𝑖

′′′𝑐𝑖�̅�𝑖𝑗 =
1

30
, 

              ∑ 𝑏𝑖
′′𝑐𝑖

3 =
1

20
,     ∑ 𝑏𝑖

′′𝑎𝑖𝑗 =
1

120
,   ∑ 𝑏𝑖

′𝑐𝑖
2 =

1

60
,    ∑ 𝑏𝑖𝑐𝑖 =

1

120
, 

order 6:  ∑ 𝑏𝑖
′′′𝑐𝑖

5 =
1

6
,   ∑ 𝑏𝑖

′′′𝑎𝑖𝑗𝑐𝑗 =
1

720
,   ∑ 𝑏𝑖

′′′�̅�𝑖𝑗𝑐𝑖
2 =

1

360
,  ∑ 𝑏𝑖

′′′𝑐𝑖
2�̅�𝑖𝑗 =

1

36
,     

              ∑ 𝑏𝑖
′′′𝑐𝑖�̅�𝑖𝑗𝑐𝑗 =

1

144
,   ∑ 𝑏𝑖

′′′𝑐𝑖𝑎𝑖𝑗 =
1

144
,  ∑ 𝑏𝑖

′′𝑐𝑖
4 =

1

30
,  ∑ 𝑏𝑖

′′�̅�𝑖𝑗𝑐𝑗 =
1

720
, 

               ∑ 𝑏𝑖
′′𝑐𝑖�̅�𝑖𝑗 =

1

180
,   ∑ 𝑏𝑖

′′�̅�𝑖𝑗 =
1

720
,   ∑ 𝑏𝑖

′𝑐𝑖
3 =

1

120
,    ∑ 𝑏𝑖𝑐𝑖

2 =
1

360
,  

               ∑ 𝑏𝑖
′�̅�𝑖𝑗 =

1

720
. 

4. Derivation Embedded EDIRKTO Methods 

To solve Equation (1) numerically, the general form of the EDIRKTO technique with 𝑚-stage is 

given. The embedded pair RK method, a current research area that is constantly enhancing existing 

programs, is then developed. The derivation of 𝑝(𝑞) pairs of implicit EDIRKTO approaches is 

used to provide minimum error estimation for step size codes. Both the order 𝑝 method 

(𝐶, 𝐴, 𝑏, 𝑏′𝑇 , 𝑏′′𝑇 , 𝑏′′′𝑇) and the order 𝑞 method (𝐶, �̅�, 𝑏1
𝑇 , 𝑏′2

𝑇 , 𝑏′′3
𝑇 , 𝑏′′′4

𝑇) constitute the foundation 

for them. The embedded pair can be launched in Butcher Tabular as follows:  

 

 

 

 

𝑐1 𝑎11             𝑎12         …         𝑎1𝑠    �̅�11      �̅�12       …        �̅�1𝑠 

𝑐2 𝑎21        𝑎22          …         𝑎2𝑠  �̅�21      �̅�22       …        �̅�2𝑠 

𝑐3 𝑎31        𝑎32          …         𝑎3𝑠 �̅�31      �̅�32       …        �̅�3𝑠 

⋮   ⋮            ⋮           ⋮           ⋮   ⋮          ⋮         ⋮           ⋮ 
𝑐𝑠 𝑎𝑠1        𝑎𝑠2          …          𝑎𝑠𝑠      �̅�𝑠1       �̅�𝑠2       …         �̅�𝑠𝑠 

  𝑏1           𝑏2        …           𝑏𝑠  

  𝑏1
′            𝑏2

′        …           𝑏𝑠
′   

  𝑏1
′′          𝑏′2

′        …          𝑏𝑠
′′  

  𝑏1
′′′        𝑏′′2

′        …         𝑏𝑠
′′′  



IHJPAS. 37 (1) 2024 

378 
 

Table 2. The Embedded Pair EDIRKTO Method 

𝐶 𝐴 �̅� 
 𝑏𝑇 

𝑏′𝑇 
𝑏′′𝑇 

 
 

 𝑏′′′𝑇  

   
 𝑏1

𝑇  

 𝑏2
′𝑇 

𝑏3
′′𝑇 

𝑏4
′′′𝑇 

 

 

    The main idea is to obtain single cost error estimation for usage in step size approach values 

before generating the embedded pair of implicit EDIRKTO approaches. The methods are 

illustrated by increasing the significant pairs and local error estimates using the values step size ℎ, 

as follows [8], 

ℎ𝑛+1 = 0.9ℎ𝑛 (
𝑇𝑜𝑙

𝐿𝑇𝐸
)

1

𝑞+1
 ,                                                                                                           (8) 

Where 𝑇𝑜𝑙 refers to the required level of accuracy, and local truncation error (LTE) is performed 

at each stage. Therefore, the step will be admitted if 𝐿𝑇𝐸 ≤  𝑇𝑜𝑙, uses the local extrapolation 

method, which indicates employing more precise calculations to boost the integration and ℎ can 

be improved by using Equation (8). If  𝐿𝑇𝐸 >  𝑇𝑜𝑙, the step will be denied and the step size ℎ will 

be cut in half. Fourth-order ODEs can be solved using the EDIRKTO technique, an embedded RK-

type method. The first pair has orders 3 and 4, whereas the second pair contains orders 4 and 5. 

These methods were built utilizing sections, which confirm that the lower-order methods produced 

the most accurate error estimates while the higher-order methods were extremely accurate. 

Therefore, doubling the step size ℎ has an impact on getting the correct results. The two derivations 

for the embedded EDIRKTO4(3) and embedded EDIRKTO5(4) methods used in this work are 

shown in Tables 3 and 4.   

The 𝐴 and 𝐶 values in EDIRKTO4(3) are calculated from the 4𝑡ℎ-order solution and then obtained 

from the three-stage embedded 3𝑟𝑑-order equation. The simultaneous solution of equations, 

followed by the solution of 𝑏1
𝑇 , 𝑏′2

𝑇 , and 𝑏′′3
𝑇 while 𝑏′′′4

𝑇 has the same values as the 4𝑡ℎ-order. The 

outcomes are as follows: 

𝑏1
𝑇 =

1

6
− 𝑏′

2
𝑇

− 𝑏′
3
𝑇

, 𝑏′
2
𝑇

= 𝑏′
2
𝑇

, 𝑏′
3
𝑇

= 𝑏′
3
𝑇

, 𝑏′′′
1
𝑇

= 0, 𝑏′′′
2
𝑇

=
1

2
, 𝑏′′′

3
𝑇

=
1

2
, 𝑏′′

1
𝑇

= −
1

2
+ 𝑏′′

3
𝑇

+

𝑏′′
3
𝑇

√3 +
1

3
√3 , 𝑏′′

2
𝑇

= −2𝑏′′
3
𝑇

+ 𝑏′′
3
𝑇

√3 −
√3

3
+ 1, 𝑏′′

3
𝑇

= 𝑏′′
3
𝑇

.   

According to [12], using the minimized command in Maple to minimize the LTE, we can calculate 

the free parameters as follows:  

𝑏′2
𝑇= -0.0142606637044632, 𝑏′3

𝑇=0.644727870311524 and 𝑏′′3
𝑇= 0.105662412164532. If the value 

is optimized in fractional form, then 𝑏′2
𝑇 = −

1

100
, 𝑏′

3
𝑇

=
6

10
, and  𝑏′′3

𝑇 =
1

10
 are chosen. 

 



IHJPAS. 37 (1) 2024 

379 
 

Table 3. The Embedded Pair EDIRKTO4 (3) Method. 

1 
 −

1

4
 −

2

1000
 

1

2
−

√3

6
     

3

10
               −

1

4
                       0        0              −

2

1000
 

1

2
+

√3

6
     

3

10
                   

3

10
                 −

1

4
        0                    

9

100
         −

2

1000
 

 
−

2

10
                   

7

10
                  

6

100
 

 

 
 

2

1000
                  

2

10
                  

2

100
 

 

 
     0                   

1

4
+

√3

12
          

1

4
−

√3

12
 

 

 
     0                         

1

2
                     

1

2
 

 

 
−

2

10
                       

7

10
                 

6

100
 

 

 
−

127

300
                −

1

100
                 

3

5
 

 

 
−

2

5
+

7√3

30
          

2

5
+

7√3

30
          

1

10
   

 

 
     0                          

1

2
                      

1

2
 

 

 

 

Moreover, from the 5𝑡ℎ-order method, the coefficients of 𝐴 and 𝐶 are computed, and a three-stage 

4𝑡ℎ-order embedded approach is then derived. The solution for 𝑏1
𝑇 , 𝑏′2

𝑇 and 𝑏′′3
𝑇 is obtained by 

simultaneously solving the equations up to order 5, while 𝑏′′′3
𝑇 has the same result as the fifth order. 

The answers are attained as 

𝑏′
1
𝑇

= −
1

12
+

2

5
𝑏′

3
𝑇

−
2√6

5
𝑏′

3
𝑇

+
√6

24
, 𝑏′

2
𝑇

= −
7

5
𝑏′

3
𝑇

+
2√6

5
𝑏′

3
𝑇

−
√6

24
+

1

4
, 𝑏′

3
𝑇

= 𝑏′
3
𝑇

 , 𝑏1
𝑇 =

1

24
− 𝑏2

𝑇 −

𝑏3
𝑇 , 𝑏2

𝑇 = 𝑏2
𝑇 , 𝑏3

𝑇 = 𝑏3
𝑇 , 𝑏′′′

1
𝑇

=
1

9
, 𝑏′′′

2
𝑇

=  
4

9
−

√6

36
, 𝑏′′′

3
𝑇

=  
4

9
+

√6

36
, 𝑏′′

1
𝑇

= 0, 𝑏′′
2
𝑇

=  
1

4
+

√6

36
 , 𝑏′′

3
𝑇

=

 
1

4
−

√6

36
  . 

We determine the free parameters according to [12] by minimizing the LTE using Maple's 

minimized command as follows: 

𝑏′3
𝑇= 0.0323023077253058, 𝑏2

𝑇= -0.148646585401452 and 𝑏3
𝑇= 0.447631755184810. If the value 

is optimized in fractional form, then 𝑏′3
𝑇 =

3

100
, 𝑏2

𝑇 = −
1

10
, and  𝑏3

𝑇 =
4

10
  are chosen. 

 

 

 

 

 



IHJPAS. 37 (1) 2024 

380 
 

Table 4. The Embedded Pair EDIRKTO5(4)  Method. 

1 
 −

3

10
 −

2

100
 

2

5
−

√6

10
     

4

10
               −

3

10
                          

279

25000
−

381√6

25000
        −

2

100
 

2

5
+

√6

10
  −

1

10
                  

5

10
                 −

3

10
 −

63

1000
+

39√6

1000
               

9

100
           −

2

100
 

  

9

10000
      

1223

6000
+

1331√6

180000
       

1223

6000
−

1331√6

180000
 

 

 
 0                  

1

12
+

√6

48
                        

1

12
−

√6

48
 

 

 
     0                

1

4
+

√6

36
                           

1

4
−

√6

36
 

 

 
    

1

9
                 

4

9
−

√6

36
                           

4

9
+

√6

36
 

 

  

  −
31

120
                           −  

1

10
                             

2

5
 

 

 
−

107

1500
+

89√6

3000
        

26

125
−

89√6

3000
                

3

100
 

 

 
     0                                   

1

4
+

√6

36
                  

1

4
−

√6

36
 

 

 
    

1

9
                                    

4

9
−

√6

36
                   

4

9
+

√6

36
 

 

   

 

5. Problems Test 

The methods presented in section 4 were tested with four different problems in this section. The 

following techniques were used to carry out the numerical experiments: 

Problem 1. [15] (Inhomogeneous linear problem) 

        𝑞(4)(𝑡) = 𝑞′(𝑡) − cos (𝑡), 𝑞(0) = −
1

2
, 𝑞′(0) =

1

2
, 𝑞′′(0) =

1

2
, 𝑞′′′(0) = −

1

2
,   

where t ∈ [0,2],  

       Exact solution is 𝑞(𝑡) =
1

2
sin(𝑡) −

1

2
cos(𝑡). 

Problem 2. [15] (Homogeneous linear problem) 

 𝑞(4)(𝑡) = 𝑞2(𝑡) + (𝑞′(𝑡))2 + sin(𝑡) − 1,   𝑞(0) = 0, 𝑞′(0) = 1, 𝑞′′(0) = 0 , 𝑞′′′(0) = −1, 

Where t ∈ [0,1],  

       The exact solution is 𝑞(𝑡) = sin (t). 

Problem 3. [15] (Homogeneous linear system) 

       𝑞1
(4)(𝑡) = −

𝑒3𝑡

4
 𝑞4

′ (𝑡), 𝑞1(0) = 1, 𝑞1
′ (0) = −1, 𝑞1

′′(0) = 1, 𝑞1
′′′(0) = −1,   

        𝑞2
(4)(𝑡) = −16 𝑒−𝑡 𝑞1

′ (𝑡), 𝑞2(0) = 1, 𝑞2
′ (0) = −2, 𝑞2

′′(0) = 4, 𝑞2
′′′(0) = −8,   

       𝑞3
(4)(𝑡) = −

81 𝑒−𝑡

2
 𝑞2

′ (𝑡), 𝑞3(0) = 1, 𝑞3
′ (0) = −3, 𝑞3

′′(0) = 9, 𝑞3
′′′(0) = −27,   



IHJPAS. 37 (1) 2024 

381 
 

      𝑞4
(4)(𝑡) = −

356 𝑒−𝑡

4
 𝑞3

′ (𝑡), 𝑞4(0) = 1, 𝑞4
′ (0) = −4, 𝑞4

′′(0) = 16, 𝑞4
′′′(0) = −64,   

where t ∈ [0,1],  

The exact solution is 𝑞1(𝑡) = 𝑒−𝑡, 𝑞2(𝑡) = 𝑒−2𝑡, 𝑞3(𝑡) = 𝑒−3𝑡, 𝑞4(𝑡) = 𝑒−4𝑡. 

Problem 4. [15] (Inhomogeneous nonlinear problem) 

 𝑞(4)(𝑡) = −
15 𝑞′(𝑡)

8 𝑞6(𝑡)
, 𝑞(0) = 1, 𝑞′(0) =

1

2
, 𝑞′′(0) = −

1

4
, 𝑞′′′(0) =

3

8
,   

Where t ∈ [0,
𝜋

4
],  

   Exact solution is 𝑞(𝑡) = √𝑡 + 1. 

 

6. Numerical Results 

The tables below show the approximation outcomes for resolving issues 1-4. In the tables, the 

following acronyms will be used: 

 Tol: Tolerance. 

 F. N: the function call's number 

 STEP: succeeded strides number 

 FSTEP: failure strides number 

 TIME: enforcement time 

 EDIRKTO4(3): the embedded 4(3) method constructed in this paper 

 EDIRKTO5(4): the embedded 5(4) method constructed in this paper 

 

Table 5. Problem 1-related numerical comparisons  

FSTEP STEP TIME F.N METHOD TOL(h) 

0 7 0.068 21 EDIRKTO4(3) 
10−2 

0 5 0.029  15 EDIRKTO5(4) 

0 19 0.081 57 EDIRKTO4(3) 
10−4 

2 14 0.036  46 EDIRKTO5(4) 

1 70 0.102 212 EDIRKTO4(3) 
10−6 

2 43 0.045 133 EDIRKTO5(4) 

 

 

 

 

 

 

 



IHJPAS. 37 (1) 2024 

382 
 

Table 6. Problem 2-related numerical comparisons  

FSTEP             STEP TIME F.N METHOD TOL(h) 

0 4 0.068 12 EDIRKTO4(3) 
10−2 

0 3 0.038 9 EDIRKTO5(4) 

1 11 0.088 35 EDIRKTO4(3) 
10−4 

0 7 0.068 21 EDIRKTO5(4) 

1 34 0.096 104 EDIRKTO4(3) 
10−6 

0 23 0.077 69 EDIRKTO5(4) 

 

Table 7. Problem 3-related numerical comparisons  

FSTEP             STEP TIME F.N METHOD TOL(h) 

2 13 0.056 43 EDIRKTO4(3) 
10−2 

2 9 0.028 31 EDIRKTO5(4) 

2 47 0.066 145 EDIRKTO4(3) 
10−4 

2 25 0.039 79 EDIRKTO5(4) 

3 205 0.105 621 EDIRKTO4(3) 
10−6 

2 76 0.087 232 EDIRKTO5(4) 

 

Table 8. Problem 4-related numerical comparisons 

FSTEP             STEP TIME F.N METHOD TOL(h) 

1 4 0.080 14 EDIRKTO4(3) 
10−2 

0 2 0.056 6 EDIRKTO5(4) 

1 10 0.089 32 EDIRKTO4(3) 
10−4 

1 6 0.057 20 EDIRKTO5(4) 

2 36 0.098 112 EDIRKTO4(3) 
10−6 

1 18 0.060 56 EDIRKTO5(4) 

 

 

 

 



IHJPAS. 37 (1) 2024 

383 
 

 

Problem 1                                                         Problem2 

  

 

Problem 3                                                       Problem 4 

Figure 1. Competence curves for methods 

7. Discussion  

     In this paper, Figure 1 illustrates the enhancements made to the Implicit Type Runge-Kutta 

Method of Diagonally Embedded Pair (EDIRKTO). This was done by plotting the decimal 

logarithm for the time curve with the highest value against the logarithm of some function call 

estimates that were taken from Tables 5-8. The current study comparison of EDIRKTO4(3) and 

EDIRKTO5(4) with the four different problems. Furthermore, Figure 1 was made using the 

corresponding numerical results from Tables 5–8, respectively. In addition, as shown in Tables 5-

8, computations of the succeeded stride number (Step) and the failure stride number (FSTEP). 

Figure 1 illustrates the numerical comparison between the EDIRKTO4(3) and EDIRKTO5(4) 

outcomes from the current study. The current study is based on the Runge-Kutta Method, which 

has been previously studied by [4, 5, 8]. However, the research at hand broadened and enhanced 

the method from explicit to implicit and from direct to diagonal.  

8. Conclusion 

   The solution of fourth ODEs using diagonally embedded implicit Runge-Kutta methods  has 

been discussed in this paper. Numerical findings demonstrate that the suggested approaches are 

much more effective in terms of the number of function evaluations while solving the general 4th-

order ODEs. 



IHJPAS. 37 (1) 2024 

384 
 

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. Alomari, A.; Anakira, N.R.; Bataineh, A.S.; Hashim, I. Approximate solution of nonlinear 

system of BVP arising in fluid flow problem. Mathematical Problems in Engineering, 2013, 

2013,1-7.  

2. Jator, S.N. Numerical integrators for fourth order initial and boundary value Problems. 

International Journal of Pure and Applied Mathematics, 2008, 47(4), 563–576. 

3. Kelesoglu, O. The solution of fourth order boundary value problem arising out of the beam-

column theory using adomian decomposition method. Mathematical Problems in Engineering, 

2014, 2014,1-6.  

4. Fawzi, F. A.; Abdullah, Z.M.; Hussein, N. K. Two Embedded Pairs for Solve Directly Third- 

    Order Ordinary Differential Equation by Using Runge-Kutta Type Method (RKTGD).  Journal 

of  Physics: Conference Series, 2021, 1879(2),1-11. 

5. Hussain, K.; Ismail, F.; Senu, N. Two embedded pairs of Runge-Kutta type methods for direct 

solution of special fourth-order ordinary differential equations. Mathematical Problems in  

Engineering, 2015, 2015,1-12. 

6. Ghawadri, N.; Senu, N.; Fawzi, F.A.; Ismail, F.; Ibrahim, Z.B. Explicit Integrator of Runge-

Kutta Type for Direct Solution of 𝑢 (4) =  𝑓 (𝑥, 𝑢, 𝑢′, 𝑢′′),  Symmetry, 2019, 11(246), 1-30. 

7. Ghawadri, N.; Senu, N.; Fawzi, F.A.; Ismail, F.; Ibrahim, Z. B. Diagonally implicit Runge–

Kutta type method for directly solving special fourth-order ordinary differential equations with 

Ill- Posed problem of a beam on elastic foundation,. Algorithms, 2018,12(10), 1-10. 

8. Jumaa, M.H.; Fawzi, F.A. The Implementations of the Embedded Diagonally Implicit Type 

Runge-Kutta Method (EDITRKM) For Special Third Order of the Ordinary Differential 

Equations, Tikrit Journal of Pure Science, 2022, 27(3), 85-91. 

9. Mechee, M . Senu, N. Ismail, F . Nikouravan, Bijan ; 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, 1-9. 

10. Demba, M.A.; Senu, N.; Ramos, H.; Watthayu, W.  Development of an Efficient Diagonally 

Implicit Runge–Kutta–Nyström 5 (4) Pair for Special Second Order IVPs, Axioms, 2022, 

11(565), 1-14. 

11. Demba, M.A.; Ramos, H.; Kumam, P.; Watthayu, W.; Senu, N.; Ahmed, I. A trigonometrically 

adapted 6 (4) explicit Runge–Kutta–Nyström pair to solve oscillating systems. Mathematical 

Methods in the Applied Sciences, 2023, 46(1), 560-678. 

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

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

press, 2018, ISBN 9781315896007. 

14. Gander, W. and Dominik G. Derivation of numerical methods using computer algebra. SIAM 

review, 1999, 41(3), 577-593. 



IHJPAS. 37 (1) 2024 

385 
 

15. Ghawadri, N.G.F. Runge-Kutta type methods for solving high-order ordinary differential 

equations, PDF thesis, UPM, 2018. 

16. 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-Haitham Journal For Pure and Applied Sciences, 2022, 35(1), 92-101. 

17. Hussain, K. A. Solving oscillation problems using optimized integrator method. Italian Journal 

of Pure and Applied Mathematics, 2022, 47, 578-587. 

18. Ismail, F.; Senu, N.; Hussain, K. A sixth-order RKFD method with four-stage for directly 

solving special fourth-order ODEs. Sains Malaysiana, 2016, 45(11), 1747-1754. 

19. Allogmany, R.; Ismail, F.; Majid, Z. A.; Ibrahim, Z. B. Implicit two-point block method for 

solving fourth-order initial value problem directly with application. Mathematical Problems in 

Engineering, 2020, 2020, 1-13. 

20. Asnor, A. I.; Yatim, S. A. M.; Ibrahim, Z. B. Numerical solution of fourth-order stiff ordinary 

differential equations by implicit high order block backward differentiation formula. In AIP 

Conference Proceedings , 2021, 2423(1). AIP Publishing. 

21. Asnor, A. I.; Yatim, S. A.; M., Ibrahim, Z. B.; Zainuddin, N. High Order Block Method for 

Third Order ODEs. Computers, Materials & Continua, 2021, 67(1),1253-1267. 

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

23. Abolarin, O. E.; Akinola, L. S.; Adeyefa, E. O.; Ogunware, B. G. Implicit hybrid block methods 

for solving second, third and fourth orders ordinary differential equations directly. Italian 

Journal of Pure and Applied Mathematics, 2022, 48, 1-21. 

24. Adeyeye, O.; Omar, Z. Implicit five-step block method with generalised equidistant points for 

solving fourth order linear and non-linear initial value problems. Ain Shams Engineering 

Journal, 2019, 10(4), 881-889. 

25. Kuboye, J.; Elusakin, O. R.; Quadri, O. F. Numerical algorithms for direct solution of fourth 

order ordinary differential equations. Journal of the Nigerian Society of Physical Sciences, 

2020, 2, 218-227. 

26. Osa, A. L.; Olaoluwa, O. E. A fifth-fourth continuous block implicit hybrid method for the 

solution of third order initial value problems in ordinary differential equations. Applied and 

Computational Mathematics, 2019, 8(3), 50-57. 

27. Adeleke, J. M.; Momoh, A. L.; Duromola, M. K. One-step Hybrid Block Method for Directly 

Solving Fifth-order Initial Value Problems of Ordinary Differential Equations. Asian Research 

Journal of Mathematics, 2022, 18(1), 53-64. 

28. Jena, S. R.; Gebremedhin, G. S. Approximate solution of a fifth order ordinary differential 

equations with block method. International Journal of Computing Science and 

Mathematics, 2020, 12(4), 413-426. 

29. Ramos, H.; Momoh, A. L. Development and implementation of a tenth-order hybrid block 

method for solving fifth-order boundary value problems. Mathematical Modelling and 

Analysis, 2021, 26(2), 267-286. 

30. Adeyeye, O.; Omar, Z. Solving third order ordinary differential equations using one-step block 

method with four equidistant generalized hybrid points. IAENG International Journal of 

Applied Mathematics, 2019, 49(2),1-9.