J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 Journal of the Nigerian Society of Physical Sciences Computational study of some three-step hybrid integrators for solution of third order ordinary differential equations Mark I. Modebeia, Olumide O. Olaiyaa,, Ignatius P. Ngwongwob, aDepartment of Mathematics Programme, National Mathematical Centre, Abuja, Nigeria bDepartment of Mathematics, Federal University of Technology, Owerri, Owerri, Nigeria. Abstract Block of hybrid methods with three off-step points based on collocation technique are presented in this work for direct approximation of solution of third-order Initial and Boundary Value Problems (IVPs and BVPs). These off-step points are formulated such that they exist only on a single step at a time. Hence, these points are shifted to three positions respectively in order to obtain three block different integrators for computational analysis. These analysis includes; order of the methods, consistency, stability and convergence, global error, number of functions evaluation and CPU time. The superiority of these methods over existing methods is established numerically on different test problems in literature. DOI:10.46481/jnsps.2021.323 Keywords: Third Order Ordinary Diferential Equation, Initial Value Problem, Boundary Value Problems, Block Method, Finite Diferences, Linear Multistep Hybrid Method Article History : Received: 29 July 2021 Received in revised form: 01 September 2021 Accepted for publication: 13 September 2021 Published: 29 November 2021 c©2021 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction Third-order problems considered in this work are of the form: y′′′ = f (x, y, y′, y′′), x ∈ [a, b] (1) with appropriate initial conditions and boundary conditions re- spetively. It is assumed that; the function f is continuous in [a, b] × R3. Also, as discussed in [1, 2], the existence and uniqueness of the solution of (1) is assumed. The problem in (1) is assumed well posed and the numerical solution is the interest in this work. Numerical methods for the solution of (1) with initial condi- tions or boundary conditions and for third order singularly per- turbed boundary value problems are numerous. Some include, Email address: gmarc.ify@gmail.com, (Olumide O. Olaiya ) Collocation method [3, 4, 5, 6], Non-polynomial splines [7], Quartic Splines [8], Adomian Decomposition Method [9], Ex- ponential Quartic Spline [10], Quartic B-spline method [11], and many others. In this work, the collocation technique is employed. This method is found to be flexible and more efficient in that since it approx- imates the solution of (1) at several intra-points, its block for- mulations consists of several linear multistep methods which are required for direct solution of (1), such that it overcomes the overlapping of pieces of solutions and it is self starting. In this work we seek the numerical solution of (1) directly with three different three-step continuous hybrid block methods with three off-step points. These methods are developed using the collocation approach. 459 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 460 2. Derivation of the Methods Here, the derivation of a continuous implicit three mid intra-step hybrid block method is described, for the solution of (1) over the integration interval [a, b], πN ≡ {a = x0 < x1 < · · · < xN−1 < xN = b} with h the constant step-size, h = xi − xi−1, i = 1, 2, . . . , N. For the solution of (1) with initial condition, the method proposed in [12], three off-step points in the interval 0 < xr < xs < xt < 1 are given such that (r, s, t) = ( 38 , 5 8 , 7 8 ). Also, an optimized two-step method with three off-step points proposed in [13] is such that r, s, v are in the interval 0 < r, s, t < 2. Here, on the interval [xn, xn+3], we consider the following subintervals; [xn, xn+1], [xn+1, xn+2] and [xn+2, xn+3] where the points (r, s, u), (r′, s′, u′) and (r′′, s′′, u′′) are the assigned off-set points in each of the subinterval respectively as xn < xn+r < xn+s < xn+u < xn+1 < xn+2 < xn+3, < xn < xn+1 < xn+r′ < xn+s′ < xn+u′ < xn+2 < xn+3 and < xn < xn+1 < xn+2 < xn+r′′ < xn+s′′ < xn+u′′ < xn+3. where (r, s, u) = ( 14 , 1 2 , 3 4 ), (r ′, s′, u′) = ( 54 , 3 2 , 7 4 ) and (r ′′, s′′, u′′) = ( 94 , 5 2 , 11 4 ). Consider the approximation w(x) of y(x) given by the polynomial y(x) ' w(x) = 9∑ i=0 ai x i . (2) and the third derivative y′′′(x) ' w′′′(x) = 9∑ i=3 ρii(i − 1)(i − 2)x i−3 (3) where ai are coefficients to be determined. 2.1. Specifications 1 To obtain the method in the interval xn < xn+r < xn+s < xn+u < xn+1 < xn+2 < xn+3, we Interpolate (2) at x = xn+i, i = 0, 1, 2 and collocating (3) at the points x = xn+ i4 , i = 0, 1, . . . , 4, 8, 12, a system of 10 equations with 10 unknown ai, i = 0, 1, . . . , 9. This system can be written in matrix form as 1 xn x2n x 3 n x 4 n x 5 n x 6 n · · · x 9 n 1 xn+1 x2n+1 x 3 n+1 x 4 n+1 x 5 n+1 x 6 n+1 · · · x 9 n+1 1 xn+2 x2n+2 x 3 n+2 x 4 n+2 x 5 n+2 x 6 n+2 · · · x 9 n+2 0 0 0 6 24xn 60x2n 120x 3 n · · · 504x 6 n 0 0 0 6 24xn+ 14 60x2 n+ 14 120x3 n+ 14 · · · 504x6 n+ 14 0 0 0 6 24xn+ 12 60x2 n+ 12 120x3 n+ 12 · · · 504x6 n+ 12 0 0 0 6 24xn+ 34 60x2 n+ 34 120x3 n+ 34 · · · 504x6 n+ 34 0 0 0 6 24xn+1 60x2n+1 120x 3 n+1 · · · 504x 6 n+1 0 0 0 6 24xn+2 60x2n+2 120x 3 n+2 · · · 504x 6 n+2 0 0 0 6 24xn+3 60x2n+3 120x 3 n+3 · · · 504x 6 n+3   a0 a1 a2 a3 a4 a5 a6 a7 a8 a9  =  yn yn+1 yn+2 h3 fn h3 fn+ 14 h3 fn+ 12 h3 fn+ 34 h3 fn+1 h3 fn+2 h3 fn+3  where y( j)n+i ' y ( j)(xn+i), fn+i ' f (xn+i, yn+i, y′n+i, y ′′ n+i), On solving the above system and obtaining the coefficients ai’s (not included here) which are then substituted into the polynomial (2), the following formula is obtained: w(x) = 2∑ i=0 αi(x)yn+i + h 3  3∑ i=0 βi(x) fn+i + 3∑ i=1 β i 4 (x) fn+ i4  (4) where the α and β are continuous coefficients (which are large expressions and are not included here, but can be easily obtained with the help of Maple software). Evaluating w(x) in (4) at the points x = xn+ 14 , xn+ 12 , xn+ 34 and xn+3 and after some simplifications, we obtain the following methods: yn+ 14 = 21yn 32 + 7yn+1 16 − 3yn+2 32 + h 3 ( 34621 fn 8847360 − 18589 f n+ 14 1013760 + 25783 f n+ 12 368640 − 46109 f n+ 34 829440 + 79271 fn+1 1474560 + 3019 fn+2 2949120 − 403 fn+3 18247680 ) yn+ 12 = 3yn 8 + 3yn+1 4 − yn+2 8 + h 3 ( 4997 fn 967680 − 5459 f n+ 14 194040 + 221 f n+ 12 2520 − 341 f n+ 34 4536 + 23099 fn+1 322560 + 3083 fn+2 2257920 − 941 fn+3 31933440 ) yn+ 32 = 5yn 32 + 15yn+1 16 − 3yn+2 32 + h 3 ( 23831 fn 6193152 − 218389 f n+ 14 9934848 + 159311 f n+ 12 2580480 − 343927 f n+ 34 5806080 + 221689 fn+1 4128768 + 36983 fn+2 36126720 − 45137 fn+3 2043740160 ) yn+3 = yn − 3yn+1 + 3yn+2 − h3 ( 149 fn 1890 − 8192 f n+ 14 24255 + 104 f n+ 12 315 + 512 f n+ 34 2835 − 251 fn+1 315 − 983 fn+2 2205 − 115 fn+3 12474 )  460 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 461 If w′(x) of (4) is evaluated at x = x i 4 , i = 0(1)4, 8, 12, following formulae for approximating the first derivatives are obtained: y′n = − 3yn 2h + 2yn+1 h − yn+2 2h + h 2 ( 139 fn 5670 − 5056 f n+ 14 72765 + 1784 f n+ 12 4725 − 12352 f n+ 34 42525 + 2161 fn+1 7560 + 181 fn+2 33075 − 443 fn+3 3742200 ) y′ n+ 14 = − 5yn 4h + 3yn+1 2h − yn+2 4h + h 2 ( 916957 fn 92897280 − 175423 f n+ 14 2661120 + 1696777 f n+ 12 9676800 − 234391 f n+ 34 1555200 + 887681 fn+1 6193152 + 422549 fn+2 154828800 − 902123 fn+3 15328051200 ) y′ n+ 12 = − yn h + yn+1 h − h 2 ( 73 fn 725760 + 577 f n+ 14 72765 + 121 f n+ 12 4725 + 337 f n+ 34 42525 + 5 fn+1 48384 − fn+2 8467200 + fn+3 119750400 ) y′ n+ 32 = − 3yn 4h + yn+1 2h + yn+2 4h − h 2 ( 958907 fn 92897280 − 1048031 f n+ 14 18627840 + 1702583 f n+ 12 9676800 − 219329 f n+ 34 1555200 + 4451939 fn+1 30965760 + 2958517 fn+2 1083801600 − 128971 fn+3 2189721600 ) y′n+1 = − yn 2h + yn+2 2h − h 2 ( 467 fn 22680 − 8768 f n+ 14 72765 + 1516 f n+ 12 4725 + 14528 f n+ 34 42525 + 533 fn+1 1890 + 1447 fn+2 264600 − 221 fn+3 1871100 ) y′n+2 = yn 2h − 2yn+1 h + 3yn+2 2h + h 2 ( 607 fn 5670 − 6592 f n+ 14 10395 + 7304 f n+ 12 4725 − 10816 f n+ 34 6075 + 7951 fn+1 7560 + 208 fn+2 4725 − 2693 fn+3 3742200 ) y′n+3 = 3yn 2h − 4yn+1 h + 5yn+2 2h − h 2 ( 22889 fn 22680 − 56384 f n+ 14 10395 + 52396 f n+ 12 4725 − 442688 f n+ 34 42525 + 169 fn+1 54 − 44249 fn+2 37800 − 105641 fn+3 1871100 )  Similarly, evaluating w′′(x) of (4), at the points x = x i 4 , i = 0(1)4, 8, 12, we obtain the formulae which approximates the second derivatives: y′′n = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 113 fn 945 + 64 f n+ 14 693 + 232 f n+ 12 315 − 1472 f n+ 34 2835 + 1411 fn+1 2520 − fn+2 90 − 61 fn+3 249480 ) y′′ n+ 14 = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 23999 fn 645120 − 59981 f n+ 14 388080 + 2921 f n+ 12 3360 − 8929 f n+ 34 15120 + 53261 fn+1 92160 + 16361 fn+2 1505280 − 4967 fn+3 21288960 ) y′′ n+ 12 = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 145 fn 3456 − 6452 f n+ 14 24255 + 73 f n+ 12 105 − 1564 f n+ 34 2835 + 22973 fn+1 40320 + 1031 fn+2 94080 − 947 fn+3 3991680 ) y′′ n+ 32 = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 77173h fn 1935360 − 96433h f n+ 14 388080 + 5561h f n+ 12 10080 − 4433h f n+ 34 6480 + 374419h fn+1 645120 + 49123h fn+2 4515840 − 14941h fn+3 63866880 ) y′′n+1 = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 107 fn 2520 − 6464 f n+ 14 24255 + 64 f n+ 12 105 − 832 f n+ 34 945 − 61 fn+1 126 + 13 fn+2 1176 − fn+3 4158 ) y′′n+2 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 374 fn 945 − 56512 f n+ 14 24255 + 584 f n+ 12 105 − 17984 f n+ 34 2835 + 1243 fn+1 360 + 391 fn+2 1470 − 739 fn+3 249480 ) y′′n+3 = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 5165 fn 1512 − 453952 f n+ 14 24255 + 1792 f n+ 12 45 − 112064 f n+ 34 2835 − 9607 fn+1 630 − 34877 fn+2 17640 − 16573 fn+3 62370 )  2.2. Specification 2 To obtain the method in the interval < xn < xn+1 < xn+r′ < xn+s′ < xn+u′ < xn+2 < xn+3, we Interpolate (2) at x = xn+i, i = 0, 1, 2 and collocating (3) at the points x = xn+ i4 , i = 0, 4(1)8, 12, a system of 10 equations with 10 unknown ai, i = 0, 1, . . . , 9 in matrix form is given as  1 xn x2n x 3 n x 4 n x 5 n x 6 n · · · x 9 n 1 xn+1 x2n+1 x 3 n+1 x 4 n+1 x 5 n+1 x 6 n+1 · · · x 9 n+1 1 xn+2 x2n+2 x 3 n+2 x 4 n+2 x 5 n+2 x 6 n+2 · · · x 9 n+2 0 0 0 6 24xn 60x2n 120x 3 n · · · 504x 6 n 0 0 0 6 24xn+1 60x2n+1 120x 3 n+1 · · · 504x 6 n+1 0 0 0 6 24xn+ 54 60x2 n+ 54 120x3 n+ 54 · · · 504x6 n+ 54 0 0 0 6 24xn+ 32 60x2 n+ 32 120x3 n+ 32 · · · 504x6 n+ 32 0 0 0 6 24xn+ 74 60x2 n+ 74 120x3 n+ 74 · · · 504x6 n+ 74 0 0 0 6 24xn+2 60x2n+2 120x 3 n+2 · · · 504x 6 n+2 0 0 0 6 24xn+3 60x2n+3 120x 3 n+3 · · · 504x 6 n+3   a0 a1 a2 a3 a4 a5 a6 a7 a8 a9  =  yn yn+1 yn+2 h3 fn h3 fn+1 h3 fn+ 54 h3 fn+ 32 h3 fn+ 74 h3 fn+2 h3 fn+3  On solving the above system and obtaining the coefficients ai’s (not included here), are then substituted into the polynomial (2), the following formula is obtained: w(x) = 2∑ i=0 αi(x)yn+i + h 3  3∑ i=0 βi(x) fn+i + 7∑ i=5 β i 4 (x) fn+ i4  (5) where the α and β are continuous coefficients (which are large expressions and are not included here, but can be easily obtained with the help of Maple software). Evaluating w(x) in (5) at the points x = xn+ 54 , xn+ 32 , xn+ 74 and xn+3 and after some simplifications, we obtain the following methods: 461 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 462 yn+ 54 = − 3yn 32 + 15yn+1 16 + 5yn+2 32 − h 3 ( 16453 fn 27095040 + 311963 fn+1 4128768 − 588473 f n+ 54 4515840 + 240025 f n+ 32 1548288 − 361883 f n+ 74 4515840 + 38035 fn+2 2064384 − 45137 fn+3 433520640 ) yn+ 32 = − yn 8 + 3yn+1 4 + 3yn+2 8 − h 3 ( 1097 fn 1354752 + 32509 fn+1 322560 − 2999 f n+ 54 17640 + 401 f n+ 32 1890 − 373 f n+ 74 3528 + 7939 fn+2 322560 − 941 fn+3 6773760 ) yn+ 74 = − 3yn 32 + 7yn+1 16 + 21yn+2 32 − h 3 ( 37607 fn 61931520 + 111511 fn+1 1474560 − 16343 f n+ 54 129024 + 180517 f n+ 32 1105920 − 9869 f n+ 74 129024 + 54527 fn+2 2949120 − 403 fn+3 3870720 ) yn+3 = yn − 3yn+1 + 3yn+2 + h3 ( 71 fn 13230 + 316 fn+1 315 − 4864 f n+ 54 2205 + 3208 f n+ 32 945 − 4864 f n+ 74 2205 + 316 fn+2 315 + 71 fn+3 13230 )  If w′(x) of (4) is evaluated at x = x i 4 , i = 0, 4(1)8, 12, the following formulae for approximating the first derivatives are obtained: y′n = − 3yn 2h + 2yn+1 h − yn+2 2h + h 2 ( 6043 fn 198450 + 13337 fn+1 7560 − 135488 f n+ 54 33075 + 2600 f n+ 32 567 − 83648 f n+ 74 33075 + 110 fn+2 189 − 2693 fn+3 793800 ) y′n+1 = − yn 2h + yn+2 2h − h 2 ( 2573 fn 793800 + 377 fn+1 945 − 23872 f n+ 54 33075 + 332 f n+ 32 405 − 14272 f n+ 74 33075 + 149 fn+2 1512 − 221 fn+3 396900 ) y′ n+ 54 = − yn 4h − yn+1 2h + 3yn+2 4h − h 2 ( 5264363 fn 3251404800 + 6257533 fn+1 30965760 − 399871 f n+ 54 1209600 + 493205 f n+ 32 1161216 − 1789537 f n+ 74 8467200 + 304673 fn+2 6193152 − 128971 fn+3 464486400 ) y′ n+ 32 = − yn+1 h + yn+2 h − h 2 ( fn 25401600 + 23 fn+1 241920 + 263 f n+ 54 33075 + 29 f n+ 32 1134 + 263 f n+ 74 33075 + 23 fn+2 241920 + fn+3 25401600 ) y′ n+ 74 = yn 4h − 3yn+1 2h + 5yn+2 4h + h 2 ( 5265037 fn 3251404800 + 6242651 fn+1 30965760 − 2879903 f n+ 54 8467200 + 2461463 f n+ 32 5806080 − 1870343 f n+ 74 8467200 + 1508483 fn+2 30965760 − 902123 fn+3 3251404800 ) y′n+2 = yn 2h − 2yn+1 h + 3yn+2 2h + h 2 ( 643 fn 198450 + 3047 fn+1 7560 − 22208 f n+ 54 33075 + 2488 f n+ 32 2835 − 12608 f n+ 74 33075 + 97 fn+2 945 − 443 fn+3 793800 ) y′n+3 = 3yn 2h − 4yn+1 h + 5yn+2 2h + h 2 ( 3697 fn 793800 + 1972 fn+1 945 − 27584 f n+ 54 4725 + 27436 f n+ 32 2835 − 244928 f n+ 74 33075 + 24713 fn+2 7560 + 2183 fn+3 56700 ) ,  Similarly, evaluating w′′(x) of (4), at the points x = x i 4 , i = 0, 4(1)8, 12, we obtain the formulae which approximates the second derivatives: y′′n = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 278 fn 1323 + 16091 fn+1 2520 − 35008 f n+ 54 2205 + 488 f n+ 32 27 − 22336 f n+ 74 2205 + 1481 fn+2 630 − 739 fn+3 52920 ) y′′n+1 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 23 fn 3528 + 13 fn+1 18 − 1216 f n+ 54 735 + 512 f n+ 32 315 − 1984 f n+ 74 2205 + 169 fn+2 840 − fn+3 882 ) y′′ n+ 54 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 2503 fn 387072 + 523829 fn+1 645120 − 5633 f n+ 54 3920 + 9313 f n+ 32 6048 − 4357 f n+ 74 5040 + 41777 fn+2 215040 − 14941 fn+3 13547520 ) y′′ n+ 32 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 5491 fn 846720 + 32443 fn+1 40320 − 580 f n+ 54 441 + 1604 f n+ 32 945 − 1964 f n+ 74 2205 + 1601 fn+2 8064 − 947 fn+3 846720 ) y′′ n+ 74 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 5843 fn 903168 + 521837 fn+1 645120 − 3155 f n+ 54 2352 + 2671 f n+ 32 1440 − 27127 f n+ 74 35280 + 41113 fn+2 215040 − 4967 fn+3 4515840 ) y′′n+2 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 43 fn 6615 + 2021 fn+1 2520 − 64 f n+ 54 49 + 1672 f n+ 32 945 − 1216 f n+ 74 2205 + 59 fn+2 210 − 61 fn+3 52920 ) y′′n+3 = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 13 fn 1512 − 2113 fn+1 630 + 5440 f n+ 54 441 − 20288 f n+ 32 945 + 5696 f n+ 74 315 − 18619 fn+2 2520 − 2851 fn+3 13230 )  2.3. Specification 3 Finally, obtaining the method in the interval xn < xn+1 < xn+2 < xn+r′′ < xn+s′′ < xn+u′′ < xn+3, we Interpolate (2) at x = xn+i, i = 0, 1, 2 and collocating (3) at the points x = xn+ i4 , i = 0, 4, 8(1)12, a system of 10 equations with 10 unknown ai, i = 0, 1, . . . , 9 in matrix form is given as 1 xn x2n x 3 n x 4 n x 5 n x 6 n · · · x 9 n 1 xn+1 x2n+1 x 3 n+1 x 4 n+1 x 5 n+1 x 6 n+1 · · · x 9 n+1 1 xn+2 x2n+2 x 3 n+2 x 4 n+2 x 5 n+2 x 6 n+2 · · · x 9 n+2 0 0 0 6 24xn 60x2n 120x 3 n · · · 504x 6 n 0 0 0 6 24xn+1 60x2n+1 120x 3 n+1 · · · 504x 6 n+1 0 0 0 6 24xn+2 60x2n+2 120x 3 n+2 · · · 504x 6 n+2 0 0 0 6 24xn+ 94 60x2 n+ 94 120x3 n+ 94 · · · 504x6 n+ 94 0 0 0 6 24xn+ 52 60x2 n+ 52 120x3 n+ 52 · · · 504x6 n+ 52 0 0 0 6 24xn+ 114 60x2 n+ 114 120x3 n+ 114 · · · 504x6 n+ 114 0 0 0 6 24xn+3 60x2n+3 120x 3 n+3 · · · 504x 6 n+3   a0 a1 a2 a3 a4 a5 a6 a7 a8 a9  =  yn yn+1 yn+2 h3 fn h3 fn+1 h3 fn+2 h3 fn+ 94 h3 fn+ 52 h3 fn+ 114 h3 fn+3  Thus, solving the above system and obtaining the coefficients ai’s, then substituted into the polynomial (2), the following formula is obtained: w(x) = 2∑ i=0 αi(x)yn+i + h 3  3∑ i=0 βi(x) fn+i + 11∑ i=9 β i 4 (x) fn+ i4  (6) 462 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 463 where the α and β are continuous coefficients (which are large expressions and are not included here, but can be easily obtained with the help of Maple software). Evaluating w(x) in (24) at the points x = xn+ 94 , xn+ 52 , xn+ 114 and xn+3 and after some simplifications, we obtain the following methods: yn+ 94 = 5yn 32 − 9yn+1 16 + 45yn+2 32 + h 3 ( 996379 fn 681246720 + 826499 fn+1 12042240 + 97453 fn+2 1376256 + 60029 f n+ 94 1935360 − 97477 f n+ 52 860160 + 247559 f n+ 114 3311616 − 33373 fn+3 2064384 ) yn+ 52 = 3yn 8 − 5yn+1 4 + 15yn+2 8 + h 3 ( 111341 fn 31933440 + 374389 fn+1 2257920 + 14657 fn+2 64512 + 169 f n+ 94 22680 − 533 f n+ 52 2520 + 6007 f n+ 114 38808 − 6721 fn+3 193536 ) yn+ 114 = 21yn 32 − 33yn+1 16 + 77yn+2 32 + h 3 ( 10073 fn 1658880 + 6018419 fn+1 20643840 + 691801 fn+2 1474560 − 10439 f n+ 94 165888 − 21131 f n+ 52 73728 + 154817 f n+ 114 645120 − 492349 fn+3 8847360 ) yn+3 = yn − 3yn+1 + 3yn+2 + h3 ( 115 fn 12474 + 983 fn+1 2205 + 251 fn+2 315 − 512 f n+ 94 2835 − 104 f n+ 52 315 + 8192 f n+ 114 24255 − 149 fn+3 1890 )  If w′(x) of (4) is evaluated at x = x i 4 , i = 0, 4(1)8, 12, following methods for approximating the first derivatives are obtained: y′n = − 3yn 2h + 2yn+1 h − yn+2 2h + h 2 ( 39883 fn 935550 + 132803 fn+1 264600 − 4087 fn+2 945 + 454208 f n+ 94 42525 − 50056 f n+ 52 4725 + 357824 f n+ 114 72765 − 20207 fn+3 22680 ) y′n+1 = − yn 2h + yn+2 2h − h 2 ( 259 fn 48600 + 11833 fn+1 66150 − 4939 fn+2 7560 + 71872 f n+ 94 42525 − 8084 f n+ 52 4725 + 5312 f n+ 114 6615 − 1661 fn+3 11340 ) y′n+2 = yn 2h − 2yn+1 h + 3yn+2 2h + h 2 ( 4423 fn 935550 + 8219 fn+1 37800 + 22 fn+2 189 + 10688 f n+ 94 42525 − 328 f n+ 52 675 + 3008 f n+ 114 10395 − 1361 fn+3 22680 ) y′ n+ 94 = 3yn 4h − 5yn+1 2h + 7yn+2 4h + h 2 ( 106886797 fn 15328051200 + 359414603 fn+1 1083801600 + 14053789 fn+2 30965760 + 60743 f n+ 94 10886400 − 4098743 f n+ 52 9676800 + 5766623 f n+ 114 18627840 − 6451643 fn+3 92897280 ) y′ n+ 52 = yn h − 3yn+1 h + 2yn+2 h + h 2 ( 1103999 fn 119750400 + 3774721 fn+1 8467200 + 192743 fn+2 241920 − 8017 f n+ 94 42525 − 1681 f n+ 52 4725 + 23999 f n+ 114 72765 − 57289 fn+3 725760 ) y′ n+ 114 = 5yn 4h − 7yn+1 2h + 9yn+2 4h + h 2 ( 25105411 fn 2189721600 + 606913043 fn+1 1083801600 + 7056257 fn+2 6193152 − 4098337 f n+ 94 10886400 − 2296823 f n+ 52 9676800 + 6636359 f n+ 114 18627840 − 8237603 fn+3 92897280 ) y′n+3 = 3yn 2h − 4yn+1 h + 5yn+2 2h + h 2 ( 51307 fn 3742200 + 6371 fn+1 9450 + 11197 fn+2 7560 − 23872 f n+ 94 42525 − 556 f n+ 52 4725 + 4544 f n+ 114 10395 − 1063 fn+3 11340 )  Similarly, evaluating w′′(x) of (4), at the points x = x i 4 , i = 0, 4(1)8, 12, we obtain the formulae which approximates the second derivatives: y′′n = yn h2 − 2yn+1 h2 + yn+2 h2 − h ( 7999 fn 31185 + 3859 fn+1 2520 − 10109 fn+2 630 + 112576 f n+ 94 2835 − 2488 f n+ 52 63 + 12736 f n+ 114 693 − 25229 fn+3 7560 ) y′′n+1 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 1013 fn 83160 + 793 fn+1 4410 − 2231 fn+2 840 + 832 f n+ 94 135 − 1856 f n+ 52 315 + 21568 f n+ 114 8085 − 299 fn+3 630 ) y′′n+2 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 8 fn 891 + 8059 fn+1 17640 + 269 fn+2 210 − 3008 f n+ 94 2835 + 88 f n+ 52 315 + 192 f n+ 114 2695 − 55 fn+3 1512 ) y′′ n+ 94 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 52169 fn 5806080 + 2062307 fn+1 4515840 + 888467 fn+2 645120 − 39223 f n+ 94 45360 + 319 f n+ 52 1440 + 3149 f n+ 114 35280 − 75403 fn+3 1935360 ) y′′ n+ 52 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 11951 fn 1330560 + 128917 fn+1 282240 + 18367 fn+2 13440 − 692 f n+ 94 945 + 23 f n+ 52 63 + 116 f n+ 114 1617 − 1487 fn+3 40320 ) y′′ n+ 114 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 573899 fn 63866880 + 2062267 fn+1 4515840 + 59125 fn+2 43008 − 4997 f n+ 94 6480 + 1087 f n+ 52 2016 + 7899 f n+ 114 43120 − 80579 fn+3 1935360 ) y′′n+3 = yn h2 − 2yn+1 h2 + yn+2 h2 + h ( 2239 fn 249480 + 403 fn+1 882 + 3419 fn+2 2520 − 1984 f n+ 94 2835 + 128 f n+ 52 315 + 10432 f n+ 114 24255 + 11 fn+3 270 )  3. Analysis of the Method 3.1. Local truncation error and order The linear difference operators associated with the formula in (4) is of the form L[z(x); h] ≡ w(x) − 2∑ j=0 α jzn+ j + h 3  3∑ j=0 βi(x)z ′′′ n+ j + 3∑ j=1 β j 4 z′′′ n+ j4  (7) Same can be formulated for (5) and (24). The Taylor series expansion of (7) around x yields L[z(x); h] = C0z(x) + C1hz ′(x) + C2h 2z′′(x) + · · · + C ph pz( p)(x) + O(h(p+1)). (8) where the Ci are constants. Suppose the first p + 3 terms vanishes, that is C0 = C1 = C2 = · · · = C p+2 = 0 and C p+3 , 0, then L[z(x); h] = C p+3h p+3y(p+3)(x) + O(hp+4) (9) 463 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 464 Here (9) is the local Truncation Error (LTE) for (7) and equivalently for (4) and p is the order. The LTE is the amount by which the exact solution of the ODE fails to satisfy the corre- sponding difference operator. The method in (4) (and similarly for (5) and (24)) is said to be consistent if p > 1 (see [14]). The LTEs of (2.1) are Cyn+1/4 = 2351 332943851520 y(10)(xn)h 10 + O(h11) (10) Cyn+1/2 = 941 20808990720 y(10)(xn)h 10 + O(h11) (11) Cyn+3/4 = 459 4110417920 y(10)(xn)h 10 + O(h11) (12) Cyn+3 = −27 1003520 y(10)(xn)h 10 + O(h11) (13) The method that approximates the first and second derivative in (2.1) and (2.1) respectively can be computed in a similar fash- ion. The LTEs of (2.2) are Cyn+5/4 = 4667875 66588770304 y(10)(xn)h 10 + O(h11) (14) Cyn+3/2 = 20439 183500800 y(10)(xn)h 10 + O(h11) (15) Cyn+7/4 = 5519899 33973862400 y(10)(xn)h 10 + O(h11) (16) Cyn+3 = 2781 5017600 y(10)(xn)h 10 + O(h11) (17) The method that approximates the first and second derivative in (2.2) and (2.2) respectively can be computed in a similar fash- ion. Finally, The LTEs of (2.3) are Cyn+9/4 = 39621879 20552089600 y(10)(xn)h 10 + O(h11) (18) Cyn+5/2 = 10239625 4161798144 y(10)(xn)h 10 + O(h11) (19) Cyn+11/4 = 5089372651 1664719257600 y(10)(xn)h 10 + O(h11) (20) Cyn+3 = 18657 5017600 + O[h]y(10)(xn)h 10 + O(h11) (21) (22) The method that approximates the first and second derivative in (2.3) and (7) respectively can be computed in a similar fashion. The methods have order p = 7. 3.2. Zero-stability and convergence A numerical method is zero-stable if the solutions remain bounded as h → 0. Following the procedure in [6, 15], to show the zero-stability the block method (2.1)-(2.1) (similarly for (2.2)-(2.2) and (2.3)-(7) respectively) may be rewritten in a form such that y(k) n+ j4 , are on the left hand side. so that the method in matrix form becomes A j0Y ( j) τ+1 = A j 1Y ( j) τ−1 + h 3(B j1 F ( j) τ+1 + B j 0 F ( j) τ−1), j = 0, 1, 2 (23) as h → 0, (8) becomes A j0Y ( j) τ+1 = A j 1Y ( j) τ−1, j = 0, 1, 2 (24) where Y ( j) τ+1 = ( y( j)1/4, y ( j) 1/2, y ( j) 3/4, y ( j) 1 , y ( j) 2 , y ( j) 3 )T , Y j µ−1 = ( y( j) −3/4, y ( j) −1/2, y ( j) −1/4, y ( j) 1 , y ( j) 2 , y ( j) −1/4 )T , A j0 is an 18 × 18 identity matrix given by A j0 =  A10 0 0 0 A20 0 0 0 A30  with A10 = A 1 0 = A 2 0 = I6×6; A j1 =  A11 0 0 0 A21 0 0 0 A31  with A11 = A 2 1 = A 3 1 =  1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0  . The characteristic polynomial of each of the matrix A11 (simi- larly for A21 and A 3 1) is given as λ 5(λ − 1) = 0. The roots of the characteristic polynomial are λ j = 0, for j = 1, . . . , 5 and λ6 = 1. Since the roots of the characteristic polynomial are all zero but one, whose modulus is 1 (see [14]), then the method ((2.1)-(2.1)) is zero stable. The method is thus convergent by the following theorem. Theorem 3.1. Henrici [16]. A linear multistep method is said to be convergent if it is consistent (with order p ≥ 1) and it is zero-stable. 3.3. Computational procedure The methods in (2.1)-(2.1) (equivalently for the methods in (2.2)-(2.2) and (2.3)-(7) respectively) are executed in single block on the interval [xn, xn+3], where n = 0, 3, . . . , N − 3, N, the number of subintervals must be divisible by 3 in order to obtain the last points on the integration interval b = xN . Thus, the methods in (2.1)-(2.1) are put in the form Z(y) = 0 such that the system is solved by the Newtons method of the form Y j+1 = Y j − (J j)−1Z j where Y = (y′0, y ′′ 0 , y1/4, y ′ 1/4, y ′′ 1/4, y1/2, y ′ 1/2, y ′′ 1/2, y3/4, y ′ 3/4, y ′′ 3/4, y1, . . . , . . . , y ′′ N ) J is the Jacobian matrix of Z. The starting values used in the Newtons method are the approximations given by the Taylor series yn+ j4 = yn + j h 4 y ′ n + 1 2 ( j h4 )2 y′′n + 1 6 ( j h4 )3 fn y′ n+ j4 = y′n + j h 4 y ′′ n + 1 2 ( j h4 )2 fn y′′ n+ j4 = y′′n + j h 4 fn (26) for j = 0(1)4, 8, 12 464 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 465 æ æ æ æ à à à à ì ì ì ì 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 10-12 10-10 10-8 10-6 CPU time M ax E rr or ì S3HI3 à S3HI2 æ S3HI1 Figure 1. Efficiency plot for S3HI1, S3HI2 and S3HI3 showing the maximum error against the CPU time. 4. Numerical test Problems and results Numerical problems are presented to show the accuracy of the three specification of methods proposed. This comparison is done based on the maximum error against the machine cost (CPU time). This is to ascertain which of the proposed methods performed favourably and in comparison to other methods in literature. The methods used are denoted as follows: • S3HI1: Proposed method described in specification 1; • S3HI2: Proposed method described in specification 2; • S3HI3: Proposed method described in specification 3. • MAE: Maximum Absolute Error obtained in [10]; • S3HBM: 3-Step Hybrid Block Method described in [6]; • OSTBM: Step Block Method described in [3]; • TS: Total nummber of sub intervals TS on [a,b]; • N: Number of steps. Problem 1.. Consider the third-order singularly perturbed bound- ary value problem discussed in [10] −�y′′′ + y = 81�2 cos(3x) + 3� sin(3t), x ∈ [0, 1] y(0) = 0, y(1) = 3� sin(3), y′(0) = 9� (27) whose exact solution is y(x) = 3� sin(3x). Figure 1 shows the efficiency curve of Maximum error against the CPU time for the methods S3HI1, S3HI2 and S3HI3. The maximum error are determined for N = 6, 12, 24 and 48 and the CPU time for each N, for S3HI1, S3HI2 and S3HI3 re- spectively. The curve reveals that S3HI2 performed better in terms of error than S3HI1 and S3HI3 but slightly costly in terms of its CPU time. Table 1 confirms the output of the efficiency curve for this problem. The S3HI2 performed better compared to other methods. æ æ æ æ à à à à ì ì ì ì 0.00 0.05 0.10 0.15 10-17 10-15 10-13 10-11 CPU time M ax E rr or ì S3HI3 à S3HI2 æ S3HI1 Figure 2. Efficiency plot for S3HI1, S3HI2 and S3HI3 showing the maximum error against the CPU time. Problem 2.. Consider the IVP discussed in [12] y′′′ = 3 sin(x), x ∈ [0, 3] y(0) = 1, y′(0) = 0, y′′(0) = 2 (28) whose exact solution is y(x) = 3 cos(x) + x 2 2 . The Figure 2 shows the efficiency curve of S3HI1, S3HI2 and S3HI3 where the maximum error obtained in each care is plot- ted against the CPU time. The maximum error are determined for N=6,12,24 and 48 and the CPU time for each N, for S3HI1, S3HI2 and S3HI3 respectively. The efficiency curve shows that S3HI3 performed better in terms of error than S3HI1 and S3HI3 which have slightly better CPU time than S3HI2. Table 2 shows errors obtain for Problem 2 with h = 0.1 us- ing the methods S3HI1, S3HI2, S3HI3 and a method of order 8 in [12]. The table agrees with the efficiency curve which in- dicate that S3HI1 and S3HI3 have the same performance while S3HI2 performed better when compared to its counterparts and the method of order 7 in [12]. Problem 3.. Consider the IVP discussed in [5] y′′′ + 4y′ = x, x ∈ [0, 1] y(0) = y′(0) = 0, y′′(0) = 1 (29) whose exact solution is y(x) = 316 (1 cos(2x)) + x2 8 . Table 3 shows the Maximum Error obtained for different values of N at the point xN = b. The results shows the superior- ity of the proposed methods when the number of steps (subin- terval) considered are less than those used in [5]. Figure 3 is the efficiency curve showing the different methods and their in- dividual performance in the absolute errors obtained against the CPU time in each case. The methods are relatively the same in terms of errors but S3HI2 exhibits more CPU time for N = 48. Problem 4.. Consider the nonlinear BVP discussed in [3]. y′′′ + 2e−3y = 4(1 + x)−3, x ∈ [0, 1] y(0) = 0, y′(0) = 1, y′(1) = 0.5 (30) 465 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 466 Table 1. Comparison of maximum Errors (ME) for Problem 1 Methods � N = 10 N = 20 N = 40 S3HI1 1/16 1.49 × 10−8 1.53 × 10−10 1.87 × 10−12 1/32 7.78 × 10−9 8.02 × 10−11 1.01 × 10−12 1/64 4.22 × 10−9 4.41 × 10−11 5.76 × 10−13 S3HI2 1/16 3.65 × 10−9 1.35 × 10−11 6.47 × 10−14 1/32 2.08 × 10−9 7.59 × 10−12 3.73 × 10−14 1/64 1.31 × 10−9 4.70 × 10−12 2.44 × 10−14 S3HI3 1/16 4.70 × 10−8 2.74 × 10−10 2.46 × 10−12 1/32 2.60 × 10−8 1.48 × 10−10 1.34 × 10−12 1/64 1.57 × 10−8 8.62 × 10−11 7.93 × 10−13 MAE in [? ] 1/16 2.32 × 10−4 6.12 × 10−5 1.52 × 10−5 1/32 9.77 × 10−5 2.59 × 10−5 6.45 × 10−6 1/64 3.78 × 10−5 1.00 × 10−6 2.50 × 10−6 Table 2. Comparison of errors for Problem 2 x Error in S3HI1 Error in S3HI2 Error in S3HI3 Error in [12] 0.1 1.2458 × 10−16 2.1548 × 10−17 5.1274 × 10−16 4.1078 × 10−15 0.2 2.2580 × 10−15 1.7852 × 10−17 4.4770 × 10−15 1.6875 × 10−14 0.3 1.2358 × 10−15 4.2878 × 10−17 1.2258 × 10−15 5.0848 × 10−14 0.4 5.2145 × 10−15 3.2154 × 10−17 1.5551 × 10−15 1.1779 × 10−13 0.5 3.0125 × 10−15 1.3358 × 10−15 1.6712 × 10−15 2.4081 × 10−13 0.6 2.1258 × 10−15 1.2015 × 10−15 7.1255 × 10−15 4.3709 × 10−13 0.7 1.1125 × 10−14 1.3287 × 10−15 1.0021 × 10−14 7.3708 × 10−13 0.8 7.1258 × 10−14 8.2148 × 10−15 4.2014 × 10−14 1.1662 × 10−12 0.9 2.2158 × 10−14 3.2358 × 10−15 2.1584 × 10−14 1.7587 × 10−12 1.0 1.0012 × 10−14 1.9985 × 10−15 1.8877 × 10−14 2.5466 × 10−12 Table 3. Comparison of maximum errors for Problem 3 S3HI1 S3HI2 S3HI3 Method of order 7 in [5] b TS ME ME ME TS ME 5 30 1.25 × 10−13 5.58 × 10−13 5.54 × 10−13 46 1.20 × 10−10 45 5.21 × 10−13 1.12 × 10−13 5.02 × 10−13 56 3.69 × 10−11 60 2.35 × 10−13 3.87 × 10−13 1.11 × 10−12 88 2.44 × 10−12 10 60 4.25 × 10−13 1.57 × 10−13 8.41 × 10−12 61 5.54 × 10−09 75 1.22 × 10−13 1.01 × 10−13 6.63 × 10−12 91 5.04 × 10−10 90 8.14 × 10−13 1.87 × 10−13 6.74 × 10−12 136 4.53 × 10−11 15 75 1.66 × 10−12 9.81 × 10−13 3.42 × 10−12 76 2.67 × 10−08 90 2.25 × 10−12 1.77 × 10−13 1.08 × 10−12 110 2.91 × 10−09 105 1.11 × 10−12 2.31 × 10−13 9.52 × 10−11 180 1.52 × 10−10 20 90 2.58 × 10−11 3.27 × 10−13 2.57 × 10−11 91 5.29 × 10−08 105 2.78 × 10−10 1.14 × 10−13 1.21 × 10−10 129 6.54 × 10−09 120 2.36 × 10−11 2.48 × 10−13 1.18 × 10−10 204 4.19 × 10−10 Table 4. Comparison of maximum Errors (ME) for Example 5 h ME in S3HI1 ME in S3HI2 ME in S3HI3 ME in [12] 1 16 5.11541 × 10 −14 2.23457 × 10−14 1.00124 × 10−15 1.03179 × 10−11 1 32 6.55230 × 10 −15 1.20048 × 10−15 5.47895 × 10−15 3.24907 × 10−13 1 64 3.75411 × 10 −16 4.78951 × 10−17 3.47785 × 10−16 1.02789 × 10−14 466 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 467 æ æ æ æ à à à à ì ì ì ì 0.02 0.04 0.06 0.08 0.10 0.12 10-13 10-12 10-11 10-10 CPU time M ax E rr or ì S3HI3 à S3HI2 æ S3HI1 Figure 3. Efficiency plot for S3HI1, S3HI2 and S3HI3 showing the maximum error against the CPU time for Problem 3 for N = 6, 12, 24, 48. æ æ æ æ à à à à ì ì ì ì ò ò ò ò ô ô ô ô 0.0 0.1 0.2 0.3 0.4 0.5 10-17 10-15 10-13 10-11 10-9 CPU time M ax E rr or ô OSTBM ò S3HBM ì S3HI3 à S3HI2 æ S3HI1 Figure 4. Efficiency plot for S3HI1, S3HI2, S3HI3, S3BHM and OS- TBM showing the maximum error against the CPU time for Problem 4 for N=6,12,24,48. whose exact solution is y(x) = ln(1 + x). Figure 4 is the efficiency curve showing the different methods and their individual performance in the absolute errors obtained against the CPU time in each case. From Figure 4, one ob- serves that, S3BHM shows good error performance relative to OSTEM (where OSTEM is an order 8 method). The proposed S3HI1, S3HI2, S3HI3 performed well in terms of error and time as such, they out performed S3BHM and OSTEM in terms of accuracy and time cost. Example 5.. Consider the BVP discussed in [12]. y′′′ + y = (x − 4) sin(x) + (1 − x) cos(x), x ∈ [0, 1] y(0) = 0, y′(0) = 1, y′(1) = −sin 1 (31) whose exact solution is y(x) = (x − 1) sin(x). Table 4 shows the maximum error obtained using different step-sizes and compared with the an order 7 method in [12]. It clearly shows that S3HI1, S3HI2 and S3HI3 are almost equiva- lent in performance but more superior to the method cited. 5. Conclusion Block of 3 points hybrid method are proposed and applied to solve third-order linear and non linear IVPs and BVPs in ordinary differential equations. The method are such that the off-step points are shifted between two step points in order to investigate which position is more efficient. All three methods possess the same characteristics, viz-a-viz, order of accuracy and number of functions evaluations. They are found to be consistent, zero stable and convergent as well. They are less ambiguous to derive. It should also be noted that this paper proposes computational comparison of 3 classes of a family of 3 steps linear multistep method. However, less emphases are placed on these methods outperforming existing methods in literature. Albeit, they performed favourably well when com- pared to other methods in literature. References [1] J. Henderson & and K. R. Prasad, ”Existence and uniqueness of solutions of three-point boundary value problems on time scales”, Journal Nonlin- ear Science 8 (2001) 1. [2] C. P. Gupta & V. Lakshmikantham, ”Existence and uniqueness theorems for a third-order three point boundary value problem”, Journal of Nonlinear Analysis: Theory and Methods in Application 15 (1991) 949. [3] R. K. Sahi, S. N. Jator & N. A. Khan, ”Continuous Fourth Derivative Method For Third-order Boundary Value Problems”, International Journal of Pure and Applied Mathematics 85 (2013) 907. [4] S. N. Jator, ”Novel Finite Difference Schemes For Third Order Boundary Value Problems” International Journal of Pure and Applied Mathematics 53 (2009) 37. [5] O. Adeyeye & Z. Omar, ”Solving Third Order Ordinary Differen- tial Equations Using One-Step Block Method with Four Equidistant Generalized Hybrid Points”, IAENG International Journal of Applied Mathematics (IJAM) 49 (2019) 1. [6] M. I. Modebei, O. O. Olaiya & A. C. Onyekonwu, ”A 3-step fourth derivatives method for numerical integrationof third order ordinary differential equations”, Int. J. Math. Ana Opt.; Theory and Applications 7 (2021) 32. [7] S. Islam, M. A. Khan, I. A. Tirmizi & E. H. Twizell, ”Non-polynomial splines approach to the solution of a system of third-order boundary-value problems”, Applied Mathematics and Computations 168 (2005) 152. [8] P. K. Pandey, ”Solving third-order Boundary Value Problems with Quartic Splines”, SpringerPlus, 5 (2016) 1. [9] Y. Q. Hasan & S. A. Alaqel, ”Application of Adomian Decomposition Method to Solving Higher Order Singular Value Problems for Ordinary Differential Equations”, Asian Journal of Probability and Statistics 9 (2020) 28. [10] A. Khan & P. Khandelwal, ”Numerical Solution of Third Order Singu- larly Perturbed Boundary Value Problems Using Exponential Quartic Spline”, Thai Journal of Mathematics 17 (2019) 663. [11] H. K. Mishra & S. Saini, ”Quartic B-Spline Method for Solving a Singular Singularly Perturbed Third-Order Boundary Value Problems”, American Journal of Numerical Analysis 3 (2015) 18. 467 Modebei et al. / J. Nig. Soc. Phys. Sci. 3 (2021) 459–468 468 [12] Ra’ft Abdelrahim, ”Numerical solution of third order boundary value problems using one-step hybrid block method”, Ain Shams Engineering Journal 10 (2019) 179. [13] B. S. H. Kashkari & S. Alqarni, ”Optimization of two-step block method with three hybrid points for solving third order initial value problems”, Journal of Nonlinear Science and Application 12 (2019) 450. [14] J. D. Lambert, ”Computational Methods in Ordinary Differential Equa- tions”, John Wiley, New York (1973). [15] M. I. Modebei, R. B. Adeniyi, S. N. Jator & H. C. Ramos, ”A block hybrid integrator for numerically solving fourth-order Initial Value Problems”, Applied Mathematics and Computation 346 (2019) 680. [16] P. Henrici, ”Discrete Variable Methods in Ordinary Differential Equa- tions”, John Wiley & Sons, New York, USA, (1962). 468