J. Nig. Soc. Phys. Sci. 5 (2023) 865 Journal of the Nigerian Society of Physical Sciences Hybrid Block Methods with Constructed Orthogonal Basis for Solution of Third-Order Ordinary Differential Equations Folake Lois Josepha, Adeyemi Sunday Olagunjub,∗, Emmanuel Oluseye Adeyefac, Adewale Adeyemi Jamesd aDepartment of Mathematical Sciences, Bingham University Karu,Nigeria bDepartment of Mathematics, Federal University of Lafia, Nigeria cDepartment of Mathematics, Federal University Oye-Ekiti, Nigeria dDepartment of Mathematics and Statistics, American University of Nigeria, Yola, Nigeria Abstract In this work, an orthogonal polynomial with a weight function w(x) = x2 + x+1 in the interval [-1, 1] was constructed and used as the basis function to develop block methods, using collocation and interpolation approach. An efficient class of continuous and discrete numerical integration schemes of implicit hybrid form for third-order problems were developed and successfully implemented. Three different problems were solved with these schemes, and they performed favourably. The investigation, using the appropriate existing theorems, shows that the methods are consistent, zero-stable, and hence, convergent. DOI:10.46481/jnsps.2023.865 Keywords: Hybrid block, Collocation, Interpolation, Third-order ODE, Integration scheme Article History : Received: 14 June 2022 Received in revised form: 27 October 2022 Accepted for publication: 27 October 2022 Published: 22 December 2022 c© 2022 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license (https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Communicated by: J. Ndam 1. Introduction There are no analytic solutions for some initial value prob- lems of third-order Ordinary Differential Equations (ODEs) of the form y′′′ = f (x, y, y′, y′′); y(a) = α, y′(a) = β, y′′(a) = γ, x ∈ (a, b) (1) On the other hand, numerical methods for solving such involves reducing the equation to systems of first order differential equa- ∗Corresponding author tel. no: +2348032515655 Email address: olagunjuas@gmail.com (Adeyemi Sunday Olagunju) tions, this technique increases the dimension of the real prob- lem, thus involves more computations, [1]. Also, linear multi- step methods (LMMs) implemented by the predictor-corrector mode have been found to be prone to error propagation and are also known to be very expensive to implement in terms of the number of function evaluations per step. The predictors often have a lower order of accuracy than the correctors, especially when all the grid and off-grid points are used for collocation and interpolation. This disadvantage led to the development of block methods from LMMs. Apart from being self-starting, a block method does not require the development of a predictor separately. The development of block methods for solving IVPs in ODEs is the central concern of this work. Many researchers 1 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 2 have used this block method to integrate IVPs, among them are [2], [3], [4], and [5]. Some researchers also constructed orthogonal polynomials using different weight functions; see [6]. These constructed polynomials were employed as basis function for the derivation of numerical schemes. They in- clude; [7] who used orthogonal polynomials with weight func- tion w(x) = x2 for interval [−1, 1]. [8] employed orthogonal polynomials with weight function w(x) = x2 in the interval [0, 1], while [9] engaged orthogonal polynomials with weight function w(x) = x for interval [0, 1]. Virtually all the a fore- mentioned contributions are for first order differential equa- tions. [10] constructed an orthogonal basis in the interval [-1, 1] with respect to the weight function w(x) = 1 + x2 , to solve the initial value problem of a third-order differential equation. [11] constructed a block method to integrate third-order ODE, [12] developed harmonic balance methods for the analysis of chaotic dynamics in nonlinear systems, while [13] formulated a block algorithm for general third-order ordinary differential equations. Several methods have been formulated using dif- ferent polynomials and targeting varying orders of ODEs (see [14], [15]). The need to further explore this area of research with the view to solving higher order differential equations effectively is the major focus of this work. To this end, orthogonal polynomials were constructed to serve as basis functions for the develop- ment of hybrid block methods for the class of problem (1). The polynomial is characterized with w(x) = x2 + x+1 in the interval [-1, 1]. 2. Formulation of the Orthogonal Polynomial We formulate here, the basis function for the numerical scheme to be developed by constructing an orthogonal polynomial Qn(x) of degree n. As established in literature, orthogonality of two polynomials is attained if their inner product varnishes in the interval within which they are defined [9]. i. e. the inner product yields:∫ b a w(x)φm(x)φn(x)d x = hnδmn for δmn = { 0,m,n 1,m=1 where w(x) is a continuous and positive weight function in the range [a,b], and the moments µ = ∫ b a w(x)xnd x, n = 0, 1, 2, ... exist. The integral 〈φm,φn〉 = ∫ b a w(x)φm(x)φn(x)d x . is the inner product of polynomials φm and φn. When it comes to orthogonality, 〈φm,φn〉 = ∫ b a w(x)φm(x)φn(x)d x = 0, f or m , n For the construction of a befitting approximating polynomial, consider the orthogonal polynomial class {φn(x)}, which is de- fined as φn(x) = n∑ r=0 C(n)r x r which requires that: φn(1) = 1 and < φm(x),φn(x) >= 0 (m , n) . Let w(x) = 1 + x + x2 and [a, b] = [−1, 1]. Then, for n = 0, we have φ0(x) = C (0) 0 and φ0(1) = 1 = C (0) 0 Hence, φ0(x) = 1 For n = 1, we have two equations φ1(1) = C (1) 0 + C (1) 1 = 1 < φ0(x),φ1(x) >= 8 3 C(1)0 + 2 3 C(1)1 = 0 These two equations yield: φ1(x) = − 1 3 + 4 3 x For n = 2, we solve C(1)0 + C (2) 1 + C (2) 2 = 1 = φ2(1) 8 3 C(1)0 + 2 3 C(2)1 + 16 5 C(2)2 = 0 =< φ0(x),φ2(x) > 6 5 C(2)1 + 8 45 C(2)2 = 0 =< φ1(x),φ2(x) > which gives the polynomial φ2(x) = − 49 66 − 10 33 x + 45 22 x2 Following this approach, we have a class of orthogonal polyno- mials given as follows: 2 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 3 Q0(x) = 1 Q1(x) = 4 3 x − 1 3 Q2(x) = 45 22 x2 − 10 33 − 49 66 Q3(x) = 1484 419 x3 − 483 838 x2 − 930 419 x + 213 838 Q4(x) = 49427 7880 x4 − 994 985 x3 − 4347 788 x2 + 2002 2995 x + 13609 23640 Q5(x) = 169587 14882 x5 − 108625 59528 x4 − 13735 1063 x3 + 50445 29764 x2 + 127805 44646 x − 5285 25512 Q6(x) = 21640593 1029296 x6 − 865683 257324 x5 − 29995065 1029296 x4 + 511695 128662 x3 + 112595805 11322256 x2 − 2700945 2830564 x − 5509685 11322256 Q7(x) = 912217735 23254947 x7 − 2332825495 372079152 x6 − 663966303 10335532 x5 + 3345316975 372079152 x4 + 1383496345 46509894 x3 − 134149785 41342128 x2 − 315107135 93019788 x + 66696665 372079152 Q8(x) = 206238477159 2794077824 x8 − 4632559789 392917194 x7 − 876515334695 6286675104 x6 + 867587721 43657466 x5 + 1026917752861 12573350208 x4 − 3774934163 392917194 x3 − 10578992655 698519456 x2 + 462916223 392917194 x + 10792421983 25146700416 Q9(x) = 7592947736959 54306345568 x9 − 4849132400307 21722582272 x8 − 4061887568739 13576586392 x7 + 2347732919259 54306345568 x6 + 5760948177621 27153172784 x5 − 2848740151257 108612691136 x4 − 750682757643 13576586392 x3 + 278681373243 54306345568 x2 + 208683420559 54306345568 x − 34811913843 217225382272 3. Numerical Block Algorithms for Third-Order Problems The analytical solution of (1) is approximated via experimental solution of the form: Y (x) = v+w−1∑ j=0 a jφ j(x) (2) where x ∈ [a, b], v and w are the number of collocation and interpolation points respectively. The function φ j(x) is the jth degree orthogonal polynomial valid in the range of integration of [a, b]. The third derivative of (2) is given by y′′′(x) = v+w−1∑ j=3 a jφ ′′′ j (x) = f (x, y, y ′, y′′) (3) To estimate the solution of the problem given in (1), interpolation must be done at least three times. As a result, equation (2) is interpolated at (xn+w) points, and equation (3) is collocated at (xn+v) points, yielding a system of equations to be solved using the Gaussian elimination method. We will use a hybrid approach to apply this concept. 3 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 4 3.1. Development of One-step Method with xn+ 14 , xn+ 12 and xn+ 34 as the Off-step Point The new off-step points are xn+ 14 , xn+ 12 and xn+ 34 . Interpolating (2)at x = xn+w, w = 0, 1 4 and 1 2 and collocating (3) at x = xn+v, v = 0, 14 , 1 2 , 3 4 and 1 ; yields the system of equations: 1 −53 53 33 −689 419 983 591 −519 310 631 375 −10287 6091 1 −1 −788 823 1059 −400 641 −159 2278 593 1002 −201 409 1 −13 −49 66 213 838 1297 2253 −301 1453 −382 785 197 1099 0 0 0 35701210 −246792 197 140399 27 −112087 7 203929 5 0 0 0 35701210 −24072 37 30440 33 8344 135 −38642 17 0 0 0 35701210 −2761 57 −17986 29 11263 59 39979 28 0 0 0 35701210 21595 39 25171 44 −30707 63 −26670 13 0 0 0 35701210 9247 8 17997 4 105317 8 64401 2   a0 a1 a2 a3 a4 a5 a6 a7  =  yn yn+ 14 yn+ 12 h3 fn h3 fn+ 14 h3 fn+ 12 h3 fn+ 34 h3 fn+1  (4) The solution of (4) gives the values: a0 = 21 20 yn + 51 20 yn+ 12 − 13 44645 yn+ 14 + 4 44645 fnh 3 + 23 2378 fn+ 12 h 3 + 23 3087 h3 fn+ 14 + 35 36069 h3 fn+ 34 + 3 64369 h3 fn+1 a1 = 35 36 yn + 89 36 yn+ 12 − 31 9 yn+ 14 + 67 928956 h3 fn + 93 8017 h3 fn+ 12 + 325 39192 h3 fn+ 14 + 89 62116 h3 fn+ 34 + 54 971405 h3 fn+1 a2 = 44 45 yn + 44 45 yn+ 12 − 88 45 yn+ 14 + 15 148586 h3 fn + 126 14065 h3 fn+ 12 + 70 11661 h3 fn+ 14 + 6 3259 h3 fn+ 34 + 6 187829 h3 fn+1 a3 = h 3 ( 3 290044 fn + 71 21932 fn+ 12 + 55 50968 fn+ 14 + 32 20807 fn+ 34 + 8 456815 fn+1 ) a4 = h 3 ( −1 96160 fn − 12 71539 fn+ 12 − 22 30633 fn+ 14 + 33 37892 fn+ 34 + 4 157477 fn+1 ) a5 = h 3 ( 1 56502 fn − 20 35299 fn+ 12 + 14 50845 fn+ 14 + 13 55005 fn+ 34 + 7 188168 fn+1 ) a6 = h 3 ( −7 233230 fn + 4 220991 fn+ 12 + 8 148167 fn+ 14 − 13 166396 fn+ 34 + 4 110967 fn+1 ) a7 = h 3 ( 5 494258 fn + 15 247129 fn+ 12 − 10 247129 fn+ 14 − 10 247129 fn+ 34 + 5 494258 fn+1 ) (5) Substituting (5) in (3) yields the continuous implicit one-step method: ȳ(x) = 1∑ j=0 α jyn+ j + α 1 4 (x)yn+ 14 + α 12 (x)yn+ 12 + α 34 (x)yn+ 34 +h  1∑ j=0 β j(x) fn+ j + β 1 4 (x) fn+ 14 + β 12 (x) fn+ 12 + β 34 (x) fn+ 34  (6) where α j(x) and β j(x) are continuous coefficients. In (6) the parameters, α j(x) and β j(x), are obtained as : α0(t) = 2t 2 + t α 1 4 (t) = −4t2 − 4t α 1 2 (t) = 2t2 + 3t + 1 β0(t) = h 3 ( t7 2520 − t6 720 + t5 180 − t4 144 − t3 19622435024537215000 − t2 23040 4 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 5 + 13t 161280 + 1 903542973820909160000 ) β 1 4 (t) = h3 ( − t7 630 + t6 720 + t5 180 − t4 144 − t3 119312808645680590000 + 199 7103 t2 + 47 6501 t − 1 34323495691077853000 ) β 1 2 (t) = h3 ( t7 420 + t6 6654769540940698600 − t5 96 + t4 1318120659808086000 + t3 48 + 21 1280 t2 + 44 12193 − 1 917808302492149500 ) β 3 4 (t) = −h3 ( t7 630 + t6 720 − t5 180 − t4 144 + t3 4590274775345433600 + 5t2 2304 + 47t 80640 − 1 13479442579811103000 ) β1(t) = h 3 ( t7 2520 + t6 1440 − t5 2880 − t4 1152 − t3 14793289125000270000 + 7t2 23040 + 7t 86843 + 1 115497606639823160000 ) (7) where t = 2x−2xn−hh . By evaluating (6) at xn+ 34 and xn+1, the main methods are obtained as: yn+ 34 = yn − 3yn+ 14 + 3yn+ 12 + h 3 ( 1 15360 fn + 21 2560 fn+ 12 − 1 3840 fn+ 34 + 1 15360 f1 ) (8) yn+1 = 3yn − 8yn+ 14 + 6yn+ 12 + h 3 ( 1 3840 fn + 43 1920 fn+ 14 + 21 640 fn+ 12 + 13 1920 fn+ 34 + 1 3820 f1 ) (9) Differentiating (6) yields the continuous coefficients: α′0(t) = 8t + 2 h α 1 4 (t) = − 16t + 8 h α′1 2 (t) = 8t + 6 h β′0(t) = h 2 ( − t6 180 + t5 120 + t4 288 − t3 144 + t2 5537761027586915300 + t 5760 − 13 86040 ) β′1 4 (t) = h2 ( − t6 45 + t5 60 + t4 18 − t3 18 − t2 2305645383593181400 + 193t 2880 + 94 6501 ) β′1 2 (t) = h2 ( t6 30 + t5 495329736356353410 − 5t4 48 + t3 164765082476010720 + t2 8 + 21t 320 + 97 13440 ) β′3 4 (t) = h2 ( t6 45 + t5 60 − t4 18 − t3 18 + t2 686884395823310340 + 5t 576 + 47 40320 ) β′1(t) = h 2 ( t6 180 + t5 120 − t4 288 − t3 144 − t2 2514196338285239800 + 7t 5760 + 13 80640 ) (10) The second derivatives of the continuous functions are given as α′′( t) = 16 h2 5 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 6 α′′1 4 (t) = − 32 h2 α′′1 2 (t) = 16 h2 β′′0 (t) = −h ( − t5 15 + t4 12 + t3 36 − t2 24 + t 1384440256896728800 + 1 ) β′′1 4 (t) = −h ( 4t5 15 − t4 6 − 4t3 9 + t2 3 + t 576411345898295360 − 193 1440 ) β′′1 2 (t) = h ( 2t5 5 + t4 40538871860771056 − t3 12000000000 + t2 29029771673218900 + t 2 21 160 ) β′′3 4 (t) = h ( − 4t5 15 + t4 6 − 4t3 9 − t2 3 + t 171721098955827550 − 5 288 ) β′′1 (t) = h ( − t5 15 − t4 12 + t3 36 + t2 24 + t 628549084571309950 − 7 2880 ) (11) The additional methods to be coupled with the main method are obtained by evaluating the first and second derivatives of (6) at xn, xn+ 14 , xn+ 12 , xn+ 34 and xn+1 respectively to give: hy′n + 6yn − 8yn+ 14 + 2yn+ 12 = h3 ( 76 19963 fn + 245 12457 fn+ 14 − 19 4480 fn+ 12 + 29 14801 fn+ 34 − 13 36149 fn+1 ) (12) hy′ n+ 14 + 2yn + 0yn+ 14 − 2yn+ 12 = h3 ( − 76 19963 fn − 62 6527 fn+ 14 − 3 8960 fn+ 12 − 1 8064 fn+ 34 − 1 32256 fn+1 ) (13) hy′ n+ 12 − 2yn + 8yn+ 14 − 6yn+ 12 = h3 ( 13 80640 fn + 94 6501 fn+ 14 + 97 13440 fn+ 12 − 47 40320 fn+ 34 + 13 80640 fn+1 ) (14) hy′ n+ 34 − 6yn + 16yn+ 14 − 10yn+ 12 = h3 ( 15 27182 fn + 209 4679 fn+ 14 + 117 1792 fn+ 12 + 58 14347 fn+ 34 + h 32256 fn+1 ) (15) hy′n+1 − 10yn + 24yn+ 14 − 14yn+ 12 = h3 ( − 11 16128 fn + 273 3596 fn+ 14 + 324 2551 fn+ 12 + 220 3527 fn+ 34 + 104 21449 fn+1 ) (16) h2y′′n − 16yn + 32yn+ 14 − 16yn+ 12 = h3 ( − 233 2880 fn − 101 480 fn+ 14 + 31 480 fn+ 12 − 41 1440 fn+ 34 + 1 192 fn+1 ) (17) h2y′′ n+ 14 − 16yn + 32yn+ 14 − 16yn+ 12 = h3 ( 1 60 fn + 1 72 fn+ 14 − 13 480 fn+ 12 + 1 120 fn+ 34 − 1 720 fn+1 ) (18) h2y′′ n+ 12 − 16yn + 32yn+ 14 − 16yn+ 12 = h3 ( − 1 2880 fn + 193 1440 fn+ 14 + 21 160 fn+ 12 − 5 288 fn+ 34 + 7 2880 fn+1 ) (19) 6 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 7 h2y′′n+1 − 16yn + 32yn+ 14 − 16yn+ 12 = h3 ( − 1 320 fn + 209 1440 fn+ 14 + 19 96 fn+ 12 + 157 480 fn+ 34 + 239 280 fn+1 ) (20) Equations (8), (9) and (12)-(20) are solved using [22] block formula defined as Aym = hBF(ym) + E(yn) + hD fn (21) A = (ai j), B = (bi j), column vectors D = (e1...er )T ,D = (d1...dr )T , ym = (yn+1...yn+r )T and f (ym) = ( fn+1, ..., fn+r )T Thus, A, B, D and E are obtained as follows: A =  3 −3 1 0 0 0 0 0 0 0 0 0 8 −6 0 1 0 0 0 0 0 0 0 0 −8 2 0 0 0 0 0 0 0 0 0 0 0 −2 0 0 1 0 0 0 0 0 0 0 8 −6 0 0 0 1 0 0 0 0 0 0 16 −10 0 0 0 0 1 0 0 0 0 0 24 −14 0 0 0 0 0 1 0 0 0 0 32 −16 0 0 0 0 0 0 0 0 0 0 32 −16 0 0 0 0 0 0 1 0 0 0 32 −16 0 0 0 0 0 0 0 1 0 0 32 −16 0 0 0 0 0 0 0 0 1 0 32 −16 0 0 0 0 0 0 0 0 0 1  , B =  29 3840 21 2560 −1 3840 1 15360 43 1920 21 640 13 1920 1 3840 245 12457 −19 4480 29 14801 −13 36149 −62 6527 −3 8960 −1 8064 1 32256 94 6501 97 13440 −47 40320 13 80640 209 4679 117 1792 58 14347 1 32256 273 3596 324 2551 220 3527 104 21449 −101 480 31 480 −41 1440 1 192 1 72 −13 480 1 120 −1 720 193 1440 21 160 −5 288 7 2880 13 120 139 480 37 360 −1 240 209 1440 19 96 157 480 239 2880  , D =  1 15360 1 3840 76 19963 −27 55121 13 80640 15 27182 11 16128 −233 2880 1 160 −1 2880 1 288 −1 320  , E =  1 0 0 3 0 0 −6 −1 0 −2 0 0 2 0 0 6 0 0 10 0 0 16 0 −1 16 0 0 16 0 0 16 0 0 16 0 0  . Substituting A, B, D and E into (21) yields the explicit schemes: yn+ 14 = yn + 1 4 y′n + h2 32 y′′n + h 3 ( 116 73583 fn + 83 50042 − 60 62633 fn+ 12 + 15 37507 fn+ 34 − 17 233341 fn+1 ) (22) yn+ 12 = yn + 1 2 y′n + h2 8 y′′n + h 3 ( 347 42269 fn + 83 5040 − 1 168 fn+ 12 + 13 5040 fn+ 34 − 19 40320 fn+1 ) (23) yn+ 34 = yn + 3 4 y′n + 9h2 32 y′′n + h 3 ( 409 70578 fn + 164 3155 − 49 7227 fn+ 12 + 45 7168 fn+ 34 − 16 14159 fn+1 ) (24) yn+1 = yn + y ′ n + h2 2 y′′n + h 3 ( 31 840 fn + 34 315 + 1 210 fn+ 12 + 2 105 fn+ 34 − 1 504 fn+1 ) (25) y′ n+ 14 = y′n + h 4 y′′n + h 2 ( 222 1393 fn + 3 128 fn+ 14 − 47 3840 fn+ 12 + 29 5760 fn+ 34 − 7 7680 fn+1 ) (26) y′ n+ 12 = y′n + h 2 y′′n + h 2 ( 53 1440 fn + 1 10 fn+ 14 − 1 48 fn+ 12 + 1 90 fn+ 34 − 1 480 fn+1 ) (27) y′ n+ 34 = y′n + 3h 4 y′′n + h 2 ( 147 2560 fn + 117 640 fn+ 14 + 27 1280 fn+ 12 + 3 128 fn+ 34 − 9 2560 fn+1 ) (28) 7 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 8 y′n+1 = y ′ n + hy ′′ n + h 2 ( 7 90 fn + 4 15 fn+ 14 + 50793 761894 fn+ 12 + 4 45 fn+ 34 + 0 ) (29) y′′ n+ 14 = y′′n + h ( 251 2880 fn + 323 1440 fn+ 14 − 11 120 fn+ 12 + 53 1440 fn+ 34 − 19 2880 fn+1 ) (30) y′′ n+ 12 = y′′n + h ( 29 360 fn + 31 90 fn+ 14 + 1 15 fn+ 12 + 1 90 fn+ 34 − 1 360 fn+1 ) (31) y′′ n+ 34 = y′′n + h ( 27 320 fn + 51 160 fn+ 14 + 9 40 fn+ 12 + 21 160 fn+ 34 − 3 320 fn+1 ) (32) y′′n+1 = y ′′ n + h ( 7 90 fn + 16 45 fn+ 14 + 2 15 fn+ 12 + 16 45 fn+ 34 7 90 fn+1 ) (33) 3.2. Development of Two-steps Method with xn+ 13 and xn+ 23 as the Off-step Point In this section, we introduce xn+ 13 and xn+ 23 as the new off-step point. Interpolating (2) at x = xn+s, s = 0, 1 3 , 2 3 and collocating (3) at x = xn+r , r = 0, 1 3 , 2 3 , 1 and 2 yields the system: 1 −53 53 33 −689 419 983 591 −519 310 631 375 −10287 6091 1 −119 73 198 236 551 −1033 1308 841 1390 −153 1999 −1835 4279 1 −79 −41 99 4586 5741 −199 1342 −397 704 852 1805 281 1441 0 0 0 8904419 −30849 197 138449 213 45884 9 0 0 0 8904419 −7449 70 13033 51 −37143 107 0 0 0 8904419 −4555 81 6221 477 12701 90 −455 44 0 0 0 8904419 −3203 529 −8993 116 3985 167 11958 67 0 0 0 8904419 13726 95 17997 32 31266 19 68426 17   a0 a1 a2 a3 a4 a5 a6 a7  =  yn yn+ 13 y 2 3 h3 fn h3 fn+ 13 h3 f n + 23 h3 fn+1 h3 fn+2  (34) Solving this system gives the values of the unknown parameters a j, j = 0(1)7 as: a0 = 157 40 yn − 48y 5 yn+ 13 + 267 40 yn+ 23 +h3  25342525 fn + 3563685 fn+1 + 3117539 fn+2 + 1871851 fn+ 13 + fn+ 232479  a1 = 31 8 yn − 10yn+ 13 + 49 8 yn+ 23 −h3 ( 241 28690 fn − 406 3287 fn+1 − 85 34601 fn+2 − 449 3939 fn+ 13 − 247 6237 fn+ 23 ) a2 = 11 5 yn − 22 5 yn+ 13 + 11 5 yn+ 23 −h3 ( 146 15675 fn − 680 5903 fn+1 − 46 16677 fn+2 + − 166 1999 fn+ 13 + 71 4953 fn+ 23 ) a3 = 264 4033 − 79 11205 fn+1 + 8 3683 fn+2 + 253 7124 fn+ 13 − 79 1611 fn+ 23 a4 = 332 15051 − 41 11313 fn+1 + 48 35893 fn+2 + 381 25870 fn+ 13 − 220 6377 fn+ 23 a5 = 37 25260 − 31 31135 fn+1 + 1 1636 fn+2 + 63 1067 fn+ 13 − 135 19337 fn+ 23 a6 = 11 135050 − 19 9761 fn+1 − 11 58979 fn+2 − 13 9535 fn+ 13 + 34 11177 fn+ 23 a7 = 47 172075 − 41 75054 fn+1 + 10 366177 fn+2 − 121 123056 fn+ 13 + 15 122039 fn+ 23 (35) Substituting (35) in (2) yields a continuous implicit two-steps method in the form: y(x) = 2∑ j=0 α j(x)yn+ j + α 1 3 (x)yn+ 13 + α 23 (x)yn+ 23 + h  1∑ j=0 β j(x) fn+ j + β 1 3 (x) fn+ 13 + β 23  (36) 8 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 9 where α j(x) and β j(x) are continuous coefficients, and the parameters α j(x) and β j(x) are given by: α0(t) = 9t2 2 + 9t 2 + 1 α 1 3 (t) = −9t2 − 12t − 3 α 2 3 (t) = 9t2 2 + 15t 2 + 3 β0(t) = fnh 3 ( 3 280 t7 − 1 1811776004788461600 t6 − 7 240 t5 − 1 48 t4 + 1 523827403997699330 t3 + 7 2160 t2 + 9 9349 t + 1 9720 ) β 1 3 (t) = fn+ 13 h 3 ( − 27 700 t7 − 9 40 t6 27 200 t5 + 9 80 t4 − 1 119422617446588190 t3 + 61 900 t2 + 601 7560 t + 49 2700 ) β 2 3 (t) = fn+ 23 h 3 ( 27 560 t7 + 9 160 t6 − 27 160 t5 − 9 32 t4 + 1 321376444486800900 t3 + 29 144 t2 + 701 6048 t + 41 2160 ) β1(t) = fn+1h 3 ( − 3 140 t7 − 3 80 t6 + 7 120 t5 + 3 16 t4 + 1 6 t3 + 11 180 t2 + 258 35179 t − 1 4860 ) β2(t) = fn+2h 3 ( 3 2800 t7 + 3 800 t6 + 11 2400 t5 + 1 480 t4 + 1 2162636848303812900 t3 − 1 5400 t2 + 1 272160 t + 1 97200 ) (37) where t = x−xn−hh . By evaluating (36) at xn+1 and xn+2 the main methods are obtained as: yn+1 = yn − 3yn+ 13 + 3yn+ 23 + h 3 ( 1 9720 fn + 49 2700 fn+ 13 + 41 2160 fn+ 23 − 1 4860 fn+1 + 1 97200 fn+2 ) (38) yn+2 = 10yn − 24yn+ 13 + 15yn+ 23 + h 3 ( − 17 486 fn + 19 54 fn+ 13 − 1 108 fn+ 23 + 205 486 fn+1 + 11 972 fn+2 ) (39) The first derivatives of continuous functions (36) are given as: α′0(t) = 9t + 9 2h α′1 3 (t) = −18t + 12 h α′2 3 (t) = 9t + 15 2h β′0(t) = h 2 ( 3 40 t6 − 1 353797023374953600 t5 − 7 48 t4 − 1 12 t3 + 1 150940880005095680 t2+ 7 1080 t + 9 9349 ) β′1 3 (t) = h2 ( − 27 100 t6 − 27 200 t5 + 27 40 t4 + 9 20 t3 − 1 36491974884133584 t2 + 61 450 t + 601 7560 ) β′2 3 (t) = h2 ( 27 80 t6 + 27 80 t5 − 27 32 t4 − 9 8 t3 + 1 90911263767416096 t2 + 29 72 t + 701 6048 ) 9 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 10 β′n+1(t) = h 2 ( − 3 20 t6 − 9 40 t5 + 7 24 t4 + 3 4 t3 + 1 2 t2 + 11 90 t + 258 35179 ) β′n+2(t) = h 2 ( 3 400 t6 + 9 400 t5 + 11 480 t4 + 1 120 t3 + 1 758390757276836610 t2 − 1 2700 t + 1 272160 ) (40) The second derivatives are given as follows: a′′0 (t) = 9 h2 a′′1 3 (t) = − 18 h2 a′′2 3 (t) = 9 h2 β′′0 (t) = h ( 9 20 t5 − 1 109969890506378530 t4 − 7 12 t3 − 1 4 t2 + 1 75470440002547824 t + 7 1080 ) β′′1 3 (t) = −h ( 81 50 t5 + 27 40 t4 − 27 10 t3 − 27 20 t2 + 1 18245987442066792 t − 61 450 ) β′′2 3 (t) = h ( 81 40 t5 + 27 16 t4 − 27 8 t3 − 27 8 t2 + 1 45455631883708040 t + 29 72 ) β′′1 (t) = h ( − 9 10 t5 − 9 80 t4 + 7 6 t3 + 9 4 t2 + t + 11 90 ) β′′2 (t) = h ( 9 200 t5 + 9 80 t4 + 11 120 t3 + 1 40 t2 + 1 379195378638418240 t − 1 2700 ) We obtain the additional methods to be be coupled with the main methods by evaluating the first and second derivatives of (36) at xn, xn+ 13 , xn+ 23 , xn+1 and xn+2 respectively to get: hy′n + 9 2 yn − 6yn+ 13 + 3 2 yn+ 23 = h3 ( 95 13608 fn + 119 3506 fn+ 13 − 17 3024 fn+ 23 + 35 19681 fn+1 − 14 328469 fn+2 ) (41) hy′ n+ 13 + 3 2 yn − 3 2 yn+ 23 = h3 ( 43 48359 fn − 127 7560 fn+ 13 − 23 30240 fn+ 23 − 1 13608 fn+1 + 1 272160 fn+2 ) (42) hy′ n+ 23 − 3 2 yn + 6yn+ 13 − 9 2 yn+ 23 = h3 ( 13 68040 fn + 186 7109 fn+ 13 181 15120 fn+ 23 − 89 68040 fn+1 + 2 1046770 fn+2 ) (43) hy′n+1 − 9 2 yn + 12yn+ 13 − 15 2 yn+ 23 = h3 ( 9 9349 fn + 601 7560 fn+ 13 + 701 6048 fn+ 23 + 258 35179 fn+1 + 193 3170 fn+2 ) (44) hy′n+2 − 27 2 yn + 30yn+ 13 − 33 2 yn+ 23 = h3 ( − 385 26240 fn + 6979 1047 fn+ 13 − 1175 1516 fn+ 23 + 652 503 fn+1 + 1 2160 fn+2 ) (45) h2y′′ n+ 13 − 9yn + 18yn+ 13 − 9yn+ 23 = h3 ( 29 3240 fn + 7 450 fn+ 13 − 11 360 fn+ 23 + 1 162 fn+1 − 1 8100 fn+2 ) (46) 10 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 11 h2y′′ n+ 23 − 9yn + 18yn+ 13 − 9yn+ 23 = h3 ( − 1 648 fn + 331 1800 fn+ 13 + 119 720 fn+ 23 − 47 3240 fn+1 + 489 32400 fn+2 ) (47) h2y′′n+1 − 9yn + 18yn+ 13 − 9yn+ 23 = h3 ( 7 1080 fn + 61 450 fn+ 13 + 29 72 fn+ 23 + 11 90 fn+1 − 1 2700 fn+2 ) (48) h2y′′n+2 − 9yn + 18yn+ 13 − 9yn+ 23 = h3 ( − 407 1080 fn + 1261 667 fn+ 13 − 1897 720 fn+ 23 + 181 72 fn+1 + 489 1786 fn+2 ) (49) Equations (38) -(39) and (42) - (50) are solved using equation (21). A, B, D and E are obtained from the equations as: A =  3 −3 1 0 0 0 0 0 0 0 0 0 24 −15 0 1 0 0 0 0 0 0 0 0 −6 32 0 0 0 0 0 0 0 0 0 0 0 −32 0 0 1 0 0 0 0 0 0 0 6 −92 0 0 0 1 0 0 0 0 0 0 12 −152 0 0 0 0 1 0 0 0 0 0 30 −332 0 0 0 0 0 1 0 0 0 0 18 −9 0 0 0 0 0 0 0 0 0 0 18 −9 0 0 0 0 0 0 1 0 0 0 18 −9 0 0 0 0 0 0 0 1 0 0 18 −9 0 0 0 0 0 0 0 0 1 0 18 −9 0 0 0 0 0 0 0 0 0 1  B =  49 2700 41 2160 −1 4860 1 97200 19 54 −1 108 205 486 11 972 119 3506 −17 3024 35 19681 −14 328469 −127 7560 −23 30240 −1 13608 1 272160 186 7109 181 15120 −89 68040 2 104677 601 7560 701 6048 258 35179 1 272160 979 1047 −1175 1516 652 503 193 3170 −97 360 47 720 −7 360 1 2160 7 450 −11 360 1 162 −1 8100 331 1800 119 720 −47 3240 7 32400 61 640 29 72 11 90 −1 2700 1261 667 −1897 720 181 72 489 1786  D =  1 9720 −17 486 95 13608 −43 48359 13 68040 9 9349 −385 2624 −119 1080 29 3240 1 648 7 1080 −407 1080  E =  1 0 0 10 0 0 −9 2 −1 0 −3 2 0 0 3 2 0 0 9 2 0 0 27 2 0 0 9 0 −1 9 0 0 9 0 0 9 0 0 9 0 0  Inputting A, B, D and E into (21) yields: yn+ 13 = yn + 1 3 y′n + 1 18 y′′n + h 3 ( 69 18185 fn + 259 70858 fn+ 13 − 53 30240 fn+ 23 + 11 22566 fn+1 − 2 173719 fn+2 ) (50) yn+ 23 = yn + 2 3 y′n + 2 9 y′′n + h 3 ( 233 11749 fn + 176 4725 fn+ 13 − 61 5670 fn+ 23 + 16 5103 fn+1 − 2 255150 fn+2 ) (51) yn+1 = yn + y ′ n + 1 2 y′′n + h 3 ( 27 560 fn + 333 2800 fn+ 13 − 9 1120 fn+ 23 + 13 1680 fn+1 − 1 5600 fn+2 ) (52) yn+2 = yn + 2y ′ n + 2y ′′ n + h 3 ( 6 35 fn + 144 175 fn+ 13 − 9 70 fn+ 23 + 16 35 fn+1 + 11 1050 fn+2 ) (53) y′ n+ 13 = +y′n + 1 3 y′′n + h 2 ( 187 6480 fn + 211 5400 fn+ 13 − 73 4320 fn+ 23 + 1 216 fn+1 − 7 64800 fn+2 ) (54) 11 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 12 y′ n+ 23 = y′n + 2 3 y′′n + h 2 ( 1 15 fn + 116 675 fn+ 13 − 7 270 fn+ 23 + 4 405 fn+1 − 1 4050 fn+2 ) (55) y′n+1 = y ′ n + y ′′ n + h 2 ( 5 48 fn + 63 200 fn+ 13 + 9 160 fn+ 23 + 1 40 fn+1 − 1 2400 fn+2 ) (56) y′ n+ 13 = y′n + 2y ′′ n + h 2 ( 1 15 fn + 36 25 fn+ 13 − 9 10 fn+ 23 + 4 3 fn+1 + 3 50 fn+2 ) (57) y′′n = y ′′ n + h ( 193 1620 fn + 57 200 fn+ 13 − 23 240 fn+ 23 + 83 3240 fn+1 − 19 32400 fn+2 ) (58) y′′ n+ 13 = y′′n + h ( 181 1620 fn + 34 75 fn+ 13 + 1 10 fn+ 23 + 2 405 fn+1 − 1 4050 fn+2 ) (59) y′′n+1 = y ′′ n + h ( 7 60 fn + 81 200 fn+ 13 + 27 80 fn+ 23 + 17 120 fn+1 − 1 1200 fn+2 ) (60) y′′n+2 = y ′′ n − h ( 4 15 fn + 54 25 fn+ 13 − 27 10 fn+ 23 + 35 15 fn+1 + 41 150 fn+2 ) (61) 4. Analysis of the Methods In this section, a detailed examination of the performance of the discussed methods is carried out. These involve basic properties of the numerical integration scheme for ODEs, which include the order, error constant, zero stability, and consistency. The main methods derived are discrete schemes belonging to the class of LMMs of the form: k∑ j=0 α jyn+ j = h 3 k∑ j=0 β j fn+ j (62) Following [16] and [18], the Local Truncation Error(LTE) associated with (63) is defined by difference operator; L[y(x) : h] = k∑ j=0 [α jy(xn + jh) − h 3β j f (xn + jh)] (63) where y(x) is an arbitrary function, continuously differentiable on [a, b]. Expanding (63) in Taylor’s series about the point x, we obtain the expression L[y(x) : h] = coy(x) + c1hy ′ (x) + ...cp+3h p+3yp+3(x) (64) where the co, c1, c2...cp...cp+3 are obtained as follows c0 = k∑ j=0 α j (65) c1 = k∑ j=1 jα j (66) c3 = 1 3! k∑ j=1 j3α j (67) cq = 1 q! [ k∑ j=1 jqα j − q(q − 1)(q − 2)(q − 3) k∑ j=1 β j j q−3] (68) In the sense of [18], equation (63) is of order p if co = c1 = c2 = c2 = ...cp = cp+1 = cp+2 = 0 and cp+3 , 0. The cp+3 , 0 is called the error constant and cp+3hp+3yp+3(xn) is the Principal Local truncation error at the point xn. equations (8),(9) and (36) all have order p= 8 with error constants C p+3 = - 1493 845571686400 , 1 1410533428 , and −8 97197 respectively, 12 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 13 4.1. Consistency of the Methods The LMM is said to be consistent if it has order p ≥ 1 and the first and second characteristic polynomials, which are defined, respectively, as ρ(z) = k∑ j=0 α jz j (69) and σ(z) = k∑ j=0 β jz j (70) where z is the principal root, satisfy the following conditions: k∑ j=0 α j = 0 (71) ρ(1) = ρ′(1) = 0 (72) and ρ′′′(1) = 3!σ(1) (Henrichi, 1962) (73) The schemes (8), (9), and (36) are of order ρ = 8 > 1 and they have been investigated to satisfy conditions (I)-(III). Hence, the schemes are consistent. 4.2. Zero Stability of the Methods The LMM is said to be Zero-stable if no root of the first characteristic polynomial ρ(R) has modulus greater than one and if and only if every root of modulus one has multiplicity not greater than the order of the differential equation. To ana- lyze the zero-stability of the methods we present equations (22) - (25) and (51)-(54) in block form as: A0ym = hB f (ym) + A ′ynhD fn where h is a fixed mesh size within a block. In line with this, the zero stability of equations (24)-(25) gives: A0 =  1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1  A′ =  0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1  B =  83 5004 −60 62633 15 3750 −17 233341 83 5040 −1 168 13 5040 −19 40320 164 3155 −49 7227 45 7168 −16 14159 34 315 1 210 2 105 −1 504  D =  0 0 0 11673583 0 0 0 34742269 0 0 0 140970578 0 0 0 31840  ρ(R) = det ( RA0 − A′ ) (74)  R 0 0 0 0 R 0 0 0 0 R 0 0 0 0 R  −  0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1  det  R 0 0 −1 0 R 0 −1 0 0 R −1 0 0 0 R − 1  e(R) = R4 − R3 = 0 R3(R − 1) = 0 R3 = 0 or R − 1 = 0 R = 0 or R = 1 Substituting A0 and A′ in equation (75) and solving for R, the values of R are obtained as 0 and 1. The block method equa- tions (24)-(33) are zero-stable, since from (75), ρ(R) = 0, sat- isfy |R j| ≤ 1, j = 1 and for those roots with |R j| = 1,the multi- plicity does not exceed 3. The zero stability of equation(51)- (54) A0 =  1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1  A′ =  0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1  B =  259 70858 −53 30240 11 22566 −2 173719 176 4725 −61 670 16 5103 −19 255150 333 2800 −9 1120 13 1680 −1 5600 144 175 −9 70 16 35 11 1050  D =  0 0 0 6918185 0 0 0 23311749 0 0 0 27560 0 0 0 635  ρ(R) = det ( RA0 − A′ ) (75)  R 0 0 0 0 R 0 0 0 0 R 0 0 0 0 R  −  0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1  det  R 0 0 −1 0 R 0 −1 0 0 R −1 0 0 0 R − 1  e(R) = R4 − R3 = 0 R3(R − 1) = 0 13 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 14 Figure 1. Region of Absolute Stability for One-Step with Off-step Points 14 , 1 2 and 34 R3 = 0 or R − 1 = 0 R = 0 or R = 1 Substituting A0 and A′ in the equation (76) and solving for R, the values of R are obtained as 0 and 1. The block method equation (51)-(54) are zero-stable, since from (76), ρ(R) = 0, satisfy |R j| ≤ 1, j = 1 and for those roots with |R j| = 1, the multiplicity does not exceed 3. 4.3. Convergence of the Methods By the theorem of Dahlguist in [17], the necessary and suf- ficient condition for an LMM to be convergent, is that, it is con- sistent and zero-stable. The methods satisfy the two conditions stated in (70)-(74), thus the methods are convergent. 4.4. Region of Absolute Stability (RAS) For the one-step with off-step Points 14 , 1 2 and 3 4 , we have y n+ 34 + 3y n+ 14 − 3y n+ 12 − yn = h3 15360 ( fn + 126 fn+ 12 − 4 f n+ 34 + fn+1 ) h(z) = 15360 ( z 3 4 − 1 + 3z 1 4 − 3z 1 2 ) z + 126z 1 2 − 4z 3 4 + 1 h(θ) = 15360 ( ei 3 4 θ − 3ei 1 2 θ + 3ei 1 4 θ ) 126ei 1 2 θ − 4ei 3 4 θ + eiθ + 1 The RAS is shown in Figure 1 For the Two-steps with Off-step Points 13 and 2 3 , we have yn+1 + 3yn+ 13 − 3y n+ 23 − yn = h3 97200 ( 10 fn + 1764 fn+ 13 + 1845 f n+ 23 − 20 fn+1 + fn+2 ) h(z) = 97200(z + 3z 1 3 − 3z 2 3 − 1) 10 + 1764z 1 3 + 1845z 2 3 − 20z + z2 h(θ) = 97200(eiθ + 3ei 1 3 θ − 3ei 2 3 θ − 1) 10 + 1764ei 1 3 θ + 1845ei 2 3 θ − 20eiθ + ei2θ Figure 2. Region of Absolute Stability for Two-steps with Off-step Points 13 and 23 The RAS is shown in the figure below 5. Application of the Methods Problem 1 Consider the nonlinear Genesio equation that governs chaotic system [12]. x′′′(t) + Ax′′(t) + Bx′(t) = x2(t) − C x(t) x(0) = 0.2, x′(0) = −0.3, x′′(0) = 0.1, t ∈ [0, 1] where A = 1.2, B = 2.29 and C = 6 are positive constants that satisfy the condition: AB < C for the existence of the solution. Problem 2 We consider the problem of thin film flow earlier studied by [11] and sourced from [19]: y′′′ = y−k, y(0) = y′(0) = y′′(0) = 1 for case k = 2, 3. We choose h = 0.1 for its solution. Problem 3 Here, the constant coefficient homogeneous problem sourced from [13] is considered. y′′′ + y′ = 0 y(0) = 0, y′(0) = 1, y′′(0) = 2 whose analytic solution is y(x) = 2(1 − cos x) + sin x is solved with step size h = 0.1. 6. Discussion of results The tables shown above present the numerical solutions of the proposed methods. Problem one, considered by [12], has 14 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 15 Table 1. Results for Problem 1 Analytical One-step with Two-steps with Solution v = 14 , 1 2 , 3 4 v = 1 3 , 2 3 0.1 0.170440346269364 0.170440347007720 0.170440345967497 0.2 0.141582173138664 0.141582172011890 0.141582152176304 0.3 0.113282963581607 0.113282961614997 0.113282616429441 0.4 0.085554524922736 0.085554528032768 0.085553513248216 0.5 0.058543682864593 0.058543687069665 0.058541676426054 0.6 0.032510877478247 0.032510880258777 0.032507616958020 0.7 0.007806854082744 0.007806857068552 0.007802314622899 0.8 -0.015152336804258 -0.015152330175052 -0.015158060443852 0.9 -0.035911645118586 -0.035911630900664 -0.035918123711900 10 -0.054004107797261 -0.054004083268490 -0.054010808999158 Table 2. Results for Problem 2 Analytical One-step with Two-steps with [11] Solution v = 14 , 1 2 , 3 4 v = 1 3 , 2 3 0.2 1.22121001337746 1.221210004612670 1.221210004446210 1.221216123 0.4 1.48883473296637 1.488834780302660 1.488834776064700 1.488884596 0.6 1.80736134919721 1.807361398646870 1.807361383887600 1.807527232 0.8 2.17981922624938 2.179819235535390 2.179819203068140 2.180203564 10 2.60827491835218 2.608274869934920 2.608274812039350 2.609006355 Table 3. Results for Problem 3 Exact One-step with Two-steps with solution v = 14 , 1 2 , 3 4 v = 1 3 , 2 3 0.1 0.109825086090778 0.109825086090808 0.109825086091021 0.2 0.238536175112581 0.238536175016773 0.238536175109910 0.3 0.384847228410130 0.384847228107163 0.384847228371941 0.4 0.547296354302880 0.547296353668227 0.547296354184755 0.5 0.724260414823453 0.724260413721916 0.724260414556519 0.6 0.913971243575675 0.913971241864159 0.913971243082311 0.7 1.114533312668710 1.114533310199220 1.114533311853540 0.8 1.323942672205190 1.323942668827970 1.323942670968180 0.9 1.540106973086150 1.540106968652950 1.540106971317480 1.0 1.760866373071620 1.760866367438980 1.760866370661570 Table 4. Error for Problem 1 One-step with Two-steps v = 14 , 1 2 , 3 4 with v = 1 3 , 2 3 0.1 7.383560000 ×10−10 3.01867000 ×10−10 0.2 1.126774000 ×10−9 2.09623600 ×10−8 0.3 1.96661000 ×10−9 3.471521660 ×10−7 0.4 3.1100320 ×10−9 1.011674520 ×10−6 0.5 4.20507200 ×10−9 2.006438539 ×10−6 0.6 2.78053000 ×10−9 3.260520227 ×10−6 0.7 2.98580800×10−9 4.539459845 ×10−6 0.8 6.62920600 ×10−9 5.723639594 ×10−6 0.9 1.42179220 ×10−8 6.478593314 ×10−6 1.0 2.45287710 ×10−8 6.701201897 ×10−6 no exact solution and, the solution of [12] is not available for comparison. However, exact numerical solutions were gener- ated via MAPLE 18 [23], and the solutions of the two proposed methods and their errors are compared in Tables 1 and 4 respec- tively. Tables 2 and 3 present the solutions to Problems 2 and 3, respectively, and in Tables 5 and 6, error comparisons with existing methods are displayed. It is obvious that the proposed methods performed favourably. 7. Conclusion Two unique hybrid algorithms have been proposed as ini- tial value problem solvers. These self-starting algorithms effec- tively and accurately handle third-order IVPs. On investigation, the schemes were found to be convergent and accurate. The ap- proaches are highly recommended for solving third-order IVPs 15 Joseph et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 865 16 Table 5. Error for Problem 2 One-step with Two steps Error in [13] v = 14 , 1 2 , 3 4 withv = 1 3 , 2 3 Order p = 4 0.2 8.7647900 ×10−9 8.9312500 ×10−9 1.070000 ×10−6 0.4 4.73362900 ×10−8 4.30983000 ×10−8 4.130000 ×10−7 0.6 4.9449660 ×10−8 3.46903000 ×10−8 8.5100000 ×10−7 0.8 9.2860100 ×10−9 2.31812400 ×10−8 1.7100000 ×10−6 1.0 4.8417260 ×10−8 1.063128300 ×10−7 3.8600000 ×10−6 Table 6. Error for Problem 3 One-step Two-steps Error in [13] v = 14 , 1 2 , 3 4 with v = 1 3 , 2 3 0.1 −3.000000 × 10−14 −1.3100000 × 10−14 1.60880000 × 10−9 0.2 9.580800000 × 10−11 1.44400000 × 10−12 1.03870000 × 10−8 0.3 3.029670000 × 10−10 2.07030000 × 10−11 2.95720000 × 10−8 0.4 6.346530000 × 10−10 6.4005000 × 10−11 2.31470000 × 10−7 0.5 1.101537000 × 10−9 1.430840000 × 10−10 4.54200000 × 10−7 0.6 1.711516000 × 10−9 2.621400000 × 10−10 1.47460000 × 10−6 0.7 2.469490000 × 10−9 4.259500000 × 10−10 2.87340000 × 10−6 0.8 3.37722000 × 10−9 6.357000000 × 10−10 4.68260000 × 10−6 0.9 4.4332000000 × 10−9 8.89120000 × 10−10 6.92170000 × 10−6 1.0 5.6326400000 × 10−9 1.184590000 × 10−9 9.59740000 × 10−6 as a numerical scheme. The development of multi-order IVP solvers will be considered in a future study. References [1] R. A. Bun & Y. D. Vasil’yev, “A numerical method for solving differential equations of any orders”, Journal of Computer Mathematical Physics 32 (1992) 3. [2] F. L. Joseph, R. B. Adeniyi & E. O. Adeyefa, “A Collocation Technique for Hybrid Block Methods with a Constructed Orthogonal basis for Sec- ond Order Ordinary Differential Equations”, Global Journal of Pure and Applied Mathematics 14 (2018) 4. [3] E. O. Adeyefa, “A Model for Solving First, Second and Third Or- der IVPs Directly”, Int. J. Appl. Comput. Math. 7 (2021) 131. https://doi.org/10.1007/s40819-021-01075-6 [4] O. E. Abolarin, J. O. Kuboye, E. O. Adeyefa & B. O. Ogunware “New efficient numerical model for solving second , third and fourth order ordi- nary differential equations directly”, Gazi University Journal of Science 33 (2020) 4. [5] J. Sunday, G. M. Kumleng, N. M. Kamoh, J. A. Kwanamu, Y. Skwame & O. Sarjiyus, “Implicit Four-Point Hybrid Block Integrator for the Sim- ulations of Stiff Models”, J. Nig. Soc. Phys. Sci. 4 (2022) 287. DOI: 10.46481/jnsps.2022.777 [6] E. O. Adeyefa, Y. Haruna, R. O. Ajewole & R. I. Ndu, “On polynomials construction”, International Journal of Mathematical Analysis 12 (2018) 6. [7] R. B. Adeniyi & H. Yahaya, “A collocation technique based on orthogonal polynomial for construction of continuous hybrid methods”, Ilorin Jour- nal of science 2 (2014) 1. https://doi.org/10.54908/iljs.2014.01.02.014. [8] R. B. Adeniyi & M. O. Bamgbola, “Formulations of Discrete and con- tinuous Hybrid method, using orthogonal polynomials as basis function”, Journal of Nigerian Association of Mathematical Physics 29 (2015) 491. [9] F. Ekundayo & R. B. Adeniyi, “A numerical integration scheme with cer- tain orthogonal polynomial in a collocation approximation technique for ODEs”, Journal of the Nigerian Association of Mathematical Physics 28 (2014) 2. [10] R. O. Folaranmi, R. B. Adeniyi & E. O. Adeyefa, “An orthogonal based block method for the solution of ODEs”, Pacific Journal of Science and Technology 17 (2017) 2. [11] E. Momoniat & F. M. Mohammed, “Symmetry reduction and numerical solution of a third-order ODE from thin ln flow”, Mathematical and Com- putational Applications, 15 (2010) 4. [12] R. Genesio & A. Tesi, “Harmonic balance methods for the analysis of chaotic dynamics in nonlinear systems”, Automatica, 28 (1992) (3) 531. http://dx.doi.org/10.1016/0005-1098(92)90177-H. [13] T. A. Anake, A. O. Adesanya, G. J. Oghonyon & M. C. Agarana, “Block algorithm for general third order ordinary differential equations”, ICAS- TOR Journal of Mathematical Sciences 7 (2013) 2. [14] E. O. Adeyefa & J. O. Kuboye, “Derivation of New Numerical Model Capable of Solving Second and Third Order Ordinary Differential Equa- tions Directly”, IAENG International Journal of Applied Mathematics 50 (2020) 2. [15] J. O. Kuboye & E. O Adeyefa, “New development numerical for- mula for solution of first and higher order ordinary differential equations”, Journal of interdisciplinary mathematics 25 (2022) 2. https://doi.org/10.1080/09720502.2021.1925453. [16] S. O. Fatunla, Numerical methods for initial value problems in ordinary differential equations, Academic press Incorporation, Harcourt Brace, Jo- vanovich Publishers, New York (1988). [17] G. Dahlquist, Some properties of linear multistep and one leg method for ordinary differential equations. Department of Computer Science, Royal Institute of Technology, Stockholm (1979). [18] J. D. Lambert, Computational methods in ordinary differ- ential equation, John Wiley and Sons, New York (1973). https://doi.org/10.1002/zamm.19740540726 [19] K. Lee, I. Fudziah & S. Norazak, “An accurate block hybrid collocation method for third order ordinary differential equations”, Journal of Ap- plied Mathematics, Hindawi Publishing Corporation 549597 (2014) 9, https://doi.org/10.1155/2014/549597 [20] G. Dahlquist, “On one-leg multistep method”, SIAM Journal on Numeri- cal Analysis 20 (1983) 1130. [21] P. Henrichi, Discrete variable methods in ODE, John Wiley and Sons, New York (1962). [22] L. F. Shampine & H. A. Watts, “Block implicit one-step methods”, Journal of Numerical Mathematics 23 (1972) 731. https://doi.org/10.1007/BF01932819 [23] A. S. Olagunju & F. L. Joseph, “Construction of Functions for Fractional Derivatives using Matlab”, Journal of Advances in Mathematics & Com- puter Science 36 (2021) 6. 16