J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 Journal of the Nigerian Society of Physical Sciences Original Research Solving Third Order Ordinary Differential Equations Directly Using Hybrid Numerical Models J. O. Kuboye, O. F. Quadri∗, O. R. Elusakin Department of Mathematics, Federal University Oye-Ekiti, Nigeria Abstract In this work, numerical methods for solving third order initial value problems of ordinary differential equations are developed. Multi-step collocation is used in deriving the methods, where power series approximate solution is employed as a basis function. Gaussian elimination approach is considered in finding the unknown variables a j, j = 0(1)8 in interpolation and collocation equations, which are substituted into the power series to give the continuous implicit schemes. The discrete schemes and its derivatives are derived by evaluating the grid and non-grid points. These schemes are arranged in a matrix form to produce the block methods. The order of the developed methods are found to be six. The numerical results proved the efficiency of the methods over the existing methods. Keywords: Interpolation, collocation, block method, multistep method and third order initial value problem Article History : Received: 21 January 2020 Received in revised form: 09 March 2020 Accepted for publication: 16 March 2020 Published: 14 May 2020 c©2020 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction In time past, to solve third order ordinary differential equa- tion, a reduction method is employed and this involves reducing the differentials to three systems of first order ordinary differ- ential equation, and the numerical approach for first order ordi- nary differential equation is used to solve the resulting system. This technique has a lot of drawbacks like wastage of time as a result of many functions to be evaluated and the computational effort required. This is extensively discussed in Ref. [1-5]. Ad- vancing in numerical methods, direct method for solving higher order initial value problems without going through the process of reduction was developed by Fatunla [2], Kuboye and Omar ∗Corresponding author tel. no: +2347068522865 Email addresses: kubbysholly2007@yahoo.com (J. O. Kuboye), qomotolaniq@gmail.com (O. F. Quadri), opeyemielusakin21@gmail.com (O. R. Elusakin) [4], Olabode [6], Olabode [7], Omar and Suleiman [8], Sagir [9], Sunday [10] amongst others. The accuracy of the meth- ods developed by these researchers are not efficient in terms of error, as this can be improved by the introduction of off grid points within the interval of integration. Hence this paper addresses the shortcoming of low accuracy found in the existing method mentioned above. 2. Methodology 2.1. Derivation of First Block Method (FBM) This section focuses on the derivation of first block method. A power series y(x) = k+4∑ j=0 a j x j (1) 69 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 70 is considered as an approximate solution to third order ordinary differential equations of the form y ′′′ = f (x, y(x), y ′ (x), y ′′ (x)), y(x0) = y0, y ′ (x0) = y ′ 0, y ′′ (x0) = y ′′ 0 . Equation (1) is differenti- ated thrice to give: y ′′′ (x) = k+4∑ j=0 j( j − 1)( j − 2)a j x j−3 (2) where step number k=4. Interpolating (1) at x = xn+i, i = 1, 2, 9 4 and collocating its highest derivative at x = xn+m, i = 0(1)4 and 94 , the resulting equations produce a system of non-linear equations of the form k+4∑ j=0 a j x j n+i = yn+i k+4∑ j=3 j( j − 1)( j − 2)a j x j−3 n+m = fn+m (3) Solving for the unknowns constants a j, j = 0(1)8 in (3) using Gaussian elimination approach and substituting into (1) to give a continuous implicit scheme of the form k−2∑ j=1 α j(t)yn+ j + α 9 4 (t)yn+ 94 = h 3 [ k∑ j=0 β j(t) fn+ j + γ 9 4 (t) fn+ 94 ] (4) where t = x−xn+k−1h α1 α2 α 9 4  =  3 5 7 5 4 5 −6 −11 −4 32 5 48 5 16 5   t0 t1 t2  (5)  β0 β1 β2 β 9 4 β3 β4  = E  t0 t1 t2 t4 t5 t6 t7 t8  (6) where E has been defined in Appendix B. Equation (4) is dif- ferentiated twice and the coefficients are as follow:  α ′ 1 α ′ 2 α ′ 9 4  =  5 87 −11 −8 48 5 32 5   t0 t1  (7)  β ′ 0 β ′ 1 β ′ 2 β ′ 9 4 β ′ 3 β ′ 4  = F  t0 t1 t3 t4 t5 t6 t7  (8) where F has been defined in Appendix B. α ′′ 1 α ′′ 2 α ′′ 9 4  =  8 5 −8 32 5  [ t0 ] (9)  β ′′ 0 β ′′ 1 β ′′ 2 β ′′ 9 4 β ′′ 3 β ′′ 4  = G  t0 t2 t3 t4 t5 t6  (10) and G =  −573440 185794560 −1892352 185794560 430080 185794560 3153920 185794560 2580480 185794560 −585111 185794560 114688 5160960 516096 5160960 215040 5160960 −860160 5160960 −774144 5160960 262947 5160960 −573440 3440640 −3268608 3440640 −3440640 3440640 6021120 3440640 7741440 3440640 496249 3440640 114688 635040 688128 635040 860160 635040 −1146880 635040 −2064384 635040 487059 635040 −573440 15482880 −3956736 15482880 −7526400 15482880 1433600 15482880 16773120 15482880 4595169 15482880 573440 144506880 4644864 144506880 13332480 144506880 16343040 144506880 7741440 144506880 −837513 144506880  The discrete schemes and its derivatives which are generated from (4) are represented in a matrix to form a block  yn+1 yn+2 yn+ 94 yn+3 yn+4  =  1 1 1 1 1  [ yn ] +  1 2 9 4 3 4  [ hy ′ n ] +  1 2 2 81 32 9 2 8  h2y ′′ n +h3  0 0 0 0 422945360 0 0 0 0 26995670 0 0 0 0 4539264373400320 0 0 0 0 8170 0 0 0 0 60562835   fn−4 fn−3 fn−2 fn−1 fn  +h3  271 1680 −19 40 8896 19845 −1003 15120 63 11760 148 105 −64 21 8192 2835 −404 945 1 30 5137263 2621440 −144020511 36700160 943569 250880 −10281087 18350080 22497669 22497669 2349 560 −243 35 1728 245 −117 112 81 980 128 15 −416 35 262144 19845 −1408 945 8 49   fn+1 fn+2 fn+ 52 fn+3 fn+4  (11) 2.2. Derivation of Second Block Method(SBM) Equation (1) is interpolated at x = xn+ j, j = 1, 2 and 5 2 . Equation(2) is collocated at x = xn+ j, j = 0(1)4 and 5 2 . The 70 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 71 same step employed in deriving the first block method (FBM) are also used and this gives the second block of the form  yn+1 yn+2 yn+ 52 yn+3 yn+4  =  1 1 1 1 1  [ yn ] +  1 2 5 2 3 4  [ hy ′ n ] +  1 2 2 25 8 9 2 8  h2y ′′ n +h3  0 0 0 0 3794032 0 0 0 0 101210 0 0 0 0 116125147456 0 0 0 0 26192240 0 0 0 0 13663   fn−4 fn−3 fn−2 fn−1 fn  +h3  115 756 −901 3360 278 945 −283 2520 59 8640 1276 945 −12 7 256 135 −76 105 83 1890 1943125 774144 −460625 172032 75125 24192 −308125 258048 225625 3096576 81 20 −4131 1120 162 35 −99 56 243 2240 7808 945 −608 105 8192 945 −128 45 40 189   fn+1 fn+2 fn+ 52 fn+3 fn+4  (12) 3. Analysis of the Method 3.1. Order of the method The method proposed by Lambert [5] is applied in finding the order of the methods (11) and (12) where y- and f-functions are expanded in Taylor series and comparing the coefficients of h, this gives the methods to be of order [6, 6, 6, 6, 6]T 3.2. Zero Stability The Continuous implicit scheme (4) is zero-stable if the roots rn = 1, 2, ..., N of the first characteristic polynomial defined by p(r) = det[A ′ − B ′ ] satisfy |rn|≤ 1 and |r|= 1 where A ′ =  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  and B ′ =  0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1  Applying the condition above the methods are found to be zero stable. 3.3. Convergence According to Henrici [3], equation(4) converges if it is zero- stable and consistent. Therefore, since the order of the methods are greater than one and also satisfies other conditions of zero stability. Hence, the new methods (11) and (12) converged. 4. Application of the method The following problems are considered in order to examine the accuracy of the new block methods. Problem 1: Consider the initial value problem(I.V.P) y ′′′ = 3 sin x, y(0) = 1, y ′ (0) = 0, y ′′ (0) = −2, h = 0.1 Theoretical solution y(x) = 3 cos x + x 2 2 − 2 Problem 2: Consider the I.V.P y ′′′ − y ′′ + y ′ − y = 0 y(0) = 1 y ′ (0) = 0 y ′′ (0) = −1 h = 0.01 with exact solution of y(x) = cos x Problem 3: Consider the initial value problem(I.V.P) y ′′′ = exp x, y(0) = 3, y ′ (0) = 1, y ′′ (0) = 5, h = 0.1 Theoretical solution y(x) = 2 + 2x2 + exp x The following acronyms are used • ES=Exact Solution • CS=Computed Solution • EIFBM=Error in first block method with s=5/2 • EISBM=Error in second block method with s=9/4 • EI0(2009)=Error in Olabode(2009) • EIO(2013)=Error in Olabode(2013) • EIS(2014)=Error in Sagir(2014) • EISU(2018)=Error in J.Sunday (2018) 5. Discussion of Results In Table 1, the results generated from the new developed block methods (EIFBM and EISBM) are compared with EIO(2009) and EIO(2013). It is clearly shown that the accuracy in terms of error reveals the efficiency of the new derived methods over EIO(2009,2013) inspite of the high step number k considered. Furthermore, the results produced from the new numerical mod- els outperform the results in EIS(2014) and EISU(2018) as evi- dent in Table 2. Additionally, the effectiveness of the new devel- oped numerical methods is demonstrated in Table 3 as it com- pared favorably with EIO(2013) in terms of error. 71 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 72 Table 1: Comparing the new method with Olabode(2009, 2013) for solving problem 1 x-value ES CS EIFBM EISBM EIO(2009) EIO(2013) 0.1 0.990012495834077020 0.990012495833478170 6.8911543E-13 5.9885430E-13 3.4077519E-11 1.65922E-10 0.2 0.960199733523725120 0.960199733519903840 4.4015902E-12 3.8212766E-12 1.2372514E-10 4.76275E-10 0.3 0.911009467376818090 0.911009467367234980 1.0999868E-11 9.5831121E-12 1.7681812E-10 6.23182E-10 0.4 0.843182982008655380 0.843182981990707400 2.0601632E-11 1.7947976E-11 4.0865533E-10 19.9134E-10 0.5 0.757747685671118280 0.757747685638492150 3.6853520E-11 3.2626124E-11 3.7111825E-10 3.28882E-10 0.6 0.656006844729034810 0.656006844668665210 6.7268413E-11 6.0369598E-11 7.0964790E-10 1.27096E-09 0.7 0.539526561853465480 0.539526561752478040 1.1150603E-10 1.0098744E-10 7.4653450E-10 4.84653E-09 0.8 0.410120128041496560 0.410120127886882570 1.6985002E-10 1.5461399E-10 1.9585035E-09 1.09585E-08 0.9 0.269829904811992980 0.269829904583073650 2.4948449E-10 2.2891933E-10 3.8880070E-09 2.0188E-08 1.0 0.120906917604419300 0.120906917269670430 3.6226498E-10 3.3474887E-10 6.3955807E-09 3.53956E-08 1.1 -0.034211635723267797 -0.034211636195096484 5.0769700E-10 4.7182869E-10 9.5232678E-09 5.66233E-08 1.2 -0.192926736569979610 -0.192926737210326740 6.8618927E-10 6.4034714E-10 1.3169979E-08 8.35700E-08 Table 2: Comparing the new method with Sagir(2014)and J.Sunday(2018) for solving problem 2 X-VALUE ES CS EIFBM EISBM EIS(2014) EISU(2018) 0.01 0.999950000416665260 0.999950000416665370 1.1102230E-16 0.0000000E+00 1.9990E-07 1.1102E-016 0.02 0.999800006666577760 0.999800006666578310 5.5511151E-16 5.5511151E-16 1.9560E-07 1.3323E-015 0.03 0.999550033748987540 0.999550033748996200 8.6597396E-15 8.7707619E-15 1.3651E-07 9.6589E-015 0.04 0.999200106660977920 0.999200106661042750 6.4837025E-14 6.4614980E-14 2.5210E-07 3.2974E-014 0.05 0.998750260394966280 0.998750260395229290 2.6301183E-14 2.6290081E-14 1.3039E-06 8.2379E-014 Table 3: Comparing the new method with OLABODE(2013) for solving problem 3 X-VALUE ES CS EIFBM EISBM EIO (2013) 0.1 3.125170918075647700 3.125170918077170500 1.5227819e-12 1.3447021e-12 9.24352E-10 0.2 3.301402758160169700 3.301402758169862000 9.6922470e-12 8.5487173e-12 18.3983E-10 0.3 3.529858807576003300 3.529858807600271000 2.4267699e-11 2.1475266e-11 24.2400E-10 0.4 3.811824697641270600 3.811824697686721800 4.5451198e-11 4.0219827e-11 53.5873E-10 0.5 4.148721270700128200 4.148721270778516200 7.8387963e-11 7.0274453e-11 7.00128E-10 0.6 4.542118800390509700 4.542118800522103200 1.3159340e-10 1.1942358e-10 3.90509E-10 0.7 4.993752707470477500 4.993752707675188400 2.0471091e-10 1.8746782e-10 6.52952E-09 0.8 5.505540928492468600 5.505540928790510200 2.9804159e-10 2.7454572e-10 2.15075E-08 0.9 6.079603111156950000 6.079603111576208400 4.1925841e-10 3.8885162e-10 3.88430E-08 1.0 6.718281828459045500 6.718281829040118500 5.8107297e-10 5.4199667e-10 6.15410E-08 72 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 73 6. Conclusion The methods are derived via multistep collocation technique which is used to develop the block methods. The methods are zero stable, consistent, converged and of order six. The ac- curacy of the methods is better when compared with existing methods in terms of error as displayed in Tables 1,2 and 3. Hence FBM and SBM are reliable and efficient numerical mod- els for solving third order initial value problems of ordinary differential equations. References [1] D. A. Awoyemi, “Class of continuous methods for general second order initial value problems in ordinary differential equations”, International Journal of Computer Mathematics 72 (2009) 29. [2] S. O. Fatunla, Numerical methods for initial value problems in ordinary differential equations, Academic press Inc. Harcourt Brace Jovanovich Publishers, New York (1988) [3] P. Henrici, “Discrete variable methods in Ordinary Differential Equa- tions”, New York, NY; Wiley 571.91 (1962) 419. [4] J. O. Kuboye & Z. Omar, “Numerical solution of third order ordinary dif- ferential equations using a seven-step block method”, International Jour- nal of Mathematical Analysis 9 (2015) 743. [5] J. D. Lambert, “Computational Method in Ordinary Differential Equa- tion”, John Wiley and Sons, Inc 571.91 (1973) 288. [6] B. T. Olabode, “A six-step scheme for the solution of fourth order ordi- nary differential equation”, The Pacific Journal of Science and Technol- ogy 10 (2009) 143. [7] B. T. Olabode, “Block multistep method for third order ordinary differen- tial equations”, FUTA Journal of Research in Sciences 2 (2013) 194. [8] Z. Omar & M. Suleiman, “Solving higher order odes directly using parallel 2-point explicit block method”, Matematika. Pengintegrasian Teknologi Dalam Sains Matematik. Universiti Sains Malaysia 21 (2005) 15. [9] A. M. Sagir, “On the Approximate Solution of Continuous Coefficients for Solving Third Order Ordinary Differential Equations”, International Journal of Mathematical, Computational, Physical, Electrical and Com- puter Engineering 8 (2018) 67. [10] J. Sunday, “On the Oscillation Criteria and Computation of Third Order Oscillatory Differential Equations”, Communications in Mathematics and Applications 9 (2018) 615. Appendix A a0 = (xn + hs) × (2h + xn) h2(s − 1) yn+1 − (xn + hs) × (h + xn) h2(s − 2) yn+2 + 2h2 + 3hxn + xn2 h2(s2 − 3s + 2) yn+s − 2h2 + 3hxn + xn2 40320h5 s ( 3h6 s6 − 51h6 s5 +331h6 s4 − 1005h6 s3 + 1363h6 s2 − 425h6 s + 1228h5 sxn −425h5 xn + 1420h4 sxn2 − 135h4 xn2 + 616h3 sxn3 + 415h3 xn3 +116h2 sxn4 + 285h2 xn4 + 8hsxn5 + 65hxn5 + 5xn6 ) fn + 2h2 + 3hxn + xn2 10080h5(s − 1) ( 3h6 s6 − 45h6 s5 +223h6 s4 − 249h6 s3 − 1193h6 s2 + 1161h6 s − 1586h5 sxn +1161h5 xn + 258h4 sxn2 − 393h4 xn2 + 406h3 sxn3 + 9h3 xn3 +102h2 sxn4 + 183h2 xn4 + 8hsxn5 + 57hxn5 + 5xn6 ) fn+1 − 2h2 + 3hxn + xn2 6720h5(s − 2) ( 3h6 s6 − 39h6 s5 +143h6 s4 + 3h6 s3 − 277h6 s2 + 31h6 s − 228h5 sxn + 31h5 xn −92h4 sxn2 + 49h4 xn2 + 252h3 sxn3 − 89h3 xn3 + 88h2 sxn4 +109h2 xn4 + 8hsxn5 + 49hxn5 + 5xn6 ) fn+2 + 2h2 + 3hxn + xn2 336h5 s(s4 − 10s3 + 35s2 − 50s + 24) ( 13h6 s5 −h6 s6 − 57h6 s4 + 83h6 s3 + 27h6 s2 − 85h6 s − 85h5 xn −27h4 xn2 + 83h3 xn3 + 57h2 xn4 + 13hxn5 + xn6 ) fn+s + 2h2 + 3hxn + xn2 10080h5(s − 3) ( 3h6 s6 − 33h6 s5 +91h6 s4 + 3h6 s3 − 173h6 s2 + 49h6 s − 158h5 sxn + 49h5 xn −50h4 sxn2 + 15h4 xn2 + 154h3 sxn3 − 47h3 xn3 + 74h2 sxn4 +63h2 xn4 + 8hsxn5 + 41hxn5 + 5xn6 ) fn+3 − 2h2 + 3hxn + xn2 40320h5(s − 4) ( 3h6 s6 − 27h6 s5 +67h6 s4 + 3h6 s3 − 125h6 s2 + 39h6 s − 116h5 sxn + 39h5 xn −36 ∗ h4 sxn2 + 9h4 xn2 + 112h3 sxn3 − 33h3 xn3 + 60h2 sxn4 +45h2 xn4 + 8hsxn5 + 33hxn5 + 5xn6 ) fn+4 a1 = − 2h + 2xn + hs h2(s − 1) yn+1 + h + 2xn + hs (h2(s − 2) yn+2 − 3h + 2xn h2(s2 − 3s + 2) yn+s + 1 40320h5 s ( 9h7 s6 − 153h7 s5 + 993h7 s4 − 3015h7 s3 +4089h7 s2 + 1181h7 s − 850h7 + 6h6 s6 xn − 102h6 s5 xn +662h6 s4 xn − 2010h6 s3 xn + 2726h6 s2 xn + 12198h6 sxn −3090h6 xn + 20160h5 sxn2 + 14000h4 sxn3 + 6720h4 xn3 +4900h3 sxn4 + 7000h3 xn4 + 840h2 sxn5 + 2940h2 xn5 + 56hsxn6 +560hxn6 + 40xn7 ) fn − 1 10080h5(s − 1) ( 9h7 s6 − 135h7 s5 + 669h7 s4 − 747h7 s3 −3579h7 s2 + 311h7 s + 2322h7 + 6h6 s6 xn − 90h6 s5 xn +446h6 s4 xn − 498h6 s3 xn − 2386h6 s2 xn − 6162h6 sxn + 5394h6 xn +6720h4 sxn3 + 3640h3 sxn4 + 3360h3 xn4 + 756h2 sxn5 +2184h2 xn5 + 56hsxn6 + 504hxn6 + 40xn7 ) fn+1 + 1 6720h5(s − 2) ( 9h7 s6 − 117h7 s5 + 429h7 s4 + 9h7 s3 −831h7 s2 − 363h7 s + 62h7 + 6h6 s6 xn − 78h6 s5 xn +286h6 s4 xn + 6h6 s3 xn − 554h6 s2 xn − 1674h6 sxn + 382h6 xn +3360h4 sxn3 + 2660h3 sxn4 + 1680h3 xn4 + 672h2 sxn5 +1596h2 xn5 + 56hsxn6 + 448hxn6 + 40xn7 ) fn+2 − 1 336h5 s(s4 − 10s3 + 35s2 − 50s + 24) ( 39h7 s5 − 3h7 s6 73 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 74 −171h7 s4 + 249h7 s3 + 81h7 s2 − 255h7 s − 170h7 − 2h6 s6 xn +26h6 s5 xn − 114h6 s4 xn + 166h6 s3 xn + 54h6 s2 xn − 170h6 sxn −618h6 xn + 1344h4 xn3 + 1400h3 xn4 + 588h2 xn5 +112hxn6 + 8xn7 ) fn+s − 1 10080h5(s − 3) ( 9h7 s6 − 99h7 s5 + 273h7 s4 + 9h7 s3 −519h7 s2 − 169h7 s + 98h7 + 6h6 s6 xn − 66h6 s5 xn +182h6 s4 xn + 6h6 s3 xn − 346h6 s2 xn − 1050h6 sxn + 354h6 xn +2240h4 sxn3 + 1960h3 sxn4 + 1120h3 xn4 + 588h2 sxn5 +1176h2 xn5 + 56hsxn6 + 392hxn6 + 40xn7 ) fn+3 + 1 40320h5(s − 4) ( 9h7 s6 − 81h7 s5 + 201h7 s4 + 9h7 s3 −375h7 s2 − 115h7 s + 78h7 + 6h6 s6 xn − 54h6 s5 xn + 134h6 s4 xn +6h6 s3 xn − 250h6 s2 xn − 762h6 sxn + 270h6 xn + 1680h4 sxn3 +1540h3 sxn4 + 840h3 xn4 + 504h2 sxn5 + 924h2 xn5 +56hsxn6 + 336hxn6 + 40xn7 ) fn+4 a2 = 1 h2(s − 1) yn+1 − 1 h2(s − 2) yn+2 + 1 h2(s2 − 3s + 2) yn+s − 1 40320h5 s ( 3h6 s6 − 51h6 s5 + 331h6 s4 − 1005h6 s3 +1363h6 s2 + 6099h6 s − 1545h6 + 20160h5 sxn + 21000h4 sxn2 +10080h4 xn2 + 9800h3 sxn3 + 14000h3 xn3 + 2100h2 sxn4 +7350h2 xn4 + 168hsxn5 + 1680hxn5 + 140xn6 ) fn + 1 10080h5(s − 1) ( 3h6 s6 − 45h6 s5 + 223h6 s4 − 249h6 s3 −1193h6 s2 − 3081h6 s + 2697h6 + 10080h4 sxn2 + 7280h3 sxn3 +6720h3 xn3 + 1890h2 sxn4 + 5460h2 xn4 + 168hsxn5 + 1512hxn5 + 140xn6 ) fn+1 − 1 6720h5(s − 2) ( 3h6 s6 − 39h6 s5 + 143h6 s4 + 3h6 s3 −277h6 s2 − 837h6 s + 191h6 + 5040h4 sxn2 + 5320h3 sxn3 +3360h3 xn3 + 1680h2 sxn4 + 3990h2 xn4 + 168hsxn5 + 1344hxn5 + 140xn6 ) fn+2 + 1 336h5 s(s4 − 10s3 + 35s2 − 50s + 24) ( 13h6 s5 − h6 s6 − 57h6 s4 + 83h6 s3 +27h6 s2 − 85h6 s − 309h6 + 2016h4 xn2 + 2800h3 xn3 + 1470h2 xn4 + 336hxn5 + 28xn6 ) fn+s + 1 10080h5(s − 3) ( 3h6 s6 − 33h6 s5 + 91h6 s4 + 3h6 s3 −173h6 s2 − 525h6 s + 177h6 + 3360h4 sxn2 + 3920h3 sxn3 +2240h3 xn3 + 1470h2 sxn4 + 2940h2 xn4 + 168hsxn5 + 1176hxn5 + 140xn6 ) fn+3 − 1 40320h5(s − 4) ( 3h6 s6 − 27h6 s5 + 67h6 s4 + 3h6 s3 −125h6 s2 − 381h6 s + 135h6 + 2520h4 sxn2 + 3080h3 sxn3 +1680h3 xn3 + 1260h2 sxn4 + 2310h2 xn4 + 168hsxn5 + 1008hxn5 + 140xn6 ) fn+4 a3 = xn + hs 144h5 s (24h4 + 50h3 xn + 35h2 xn2 + 10hxn3 + xn4) fn 74 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 75 − xn(xn + hs) 36h5(s − 1) (24h3 + 26h2 xn + 9hxn2 + xn3) fn+1 + xn(xn + hs) 24h5(s − 2) (12h3 + 19h2 xn + 8hxn2 + xn3) fn+2 − xn 6h5 s(s4 − 10s3 + 35s2 − 50s + 24) (24h4 + 50h3 xn + 35h2 xn2 + 10hxn3 + xn4) fn+s − xn(xn + hs) 36h5(s − 3) (8h3 + 14h2 xn + 7hxn2 + xn3) fn+3 + xn(xn + hs) 144h5(s − 4) (6h3 + 11h2 xn + 6hxn2 + xn3) fn+4 a4 = − 1 576h5 s (50h4 s + 40hxn3 + 100h3 xn + 24h4 + 5xn4 + 105h2 xn2 + 4hsxn3 + 70h3 sxn + 30h2 sxn2) fn + 1 144h5(s − 1) (24h4 s + 36hxn3 + 48h3 xn + 5xn4 + 78h2 xn2 + 4hsxn3 + 52h3 sxn + 27h2 sxn2) fn+1 − 1 96h5(s − 2) (12h4 s + 32hxn3 + 24h3 xn + 5xn4 + 57h2 xn2 + 4hsxn3 + 38h3 sxn + 24h2 sxn2) fn+2 + 1 24h5 s(s4 − 10s3 + 35s2 − 50s + 24) (24h4 + 100h3 xn + 105h2 xn2 + 40hxn3 + 5xn4) fn+s + 1 144h5(s − 3) (8h4 s + 28hxn3 + 16h3 xn + 5xn4 + 42h2 xn2 + 4hsxn3 + 28h3 sxn + 21h2 sxn2)) fn+3 − 1 576h5(s − 4) (6h4 s + 24hxn3 + 12h3 xn + 5xn4 + 33h2 xn2 + 4hsxn3 + 22h3 sxn + 18h2 sxn2) fn+4 a5 = 1 1440h5 s (35h3 s + 60hxn2 + 105h2 xn + 50h3 + 10xn3 + 6hsxn2 + 30h2 sxn) fn − 1 360 ∗ h5 ∗ (s − 1) (26h3 s + 54hxn2 + 78h2 xn + 24h3 + 10xn3 + 6hsxn2 + 27h2 sxn) f n + 1 + 1 240h5(s − 2) (19h3 s + 48hxn2 + 57h2 xn + 12h3 + 10xn3 + 6hsxn2 + 24h2 sxn)) fn+2 − 1 12h5 s(s4 − 10s3 + 35s2 − 50s + 24) (10h3 + 21h2 xn + 12hxn2 + 2xn3) fn+s − 1 360h5(s − 3) (14h3 s + 42hxn2 + 42h2 xn + 8h3 + 10xn3 + 6hsxn2 + 21h2 sxn) fn+3 + 1 1440h5(s − 4) (11h3 s + 36hxn2 + 33h2 xn + 6h3 + 10xn3 + 6hsxn2 + 18h2 sxn) fn+4 a6 = − 1 2880h5 s (40hxn + 10h2 s + 35h2 + 10xn2 + 4hsxn) fn + 1 720h5(s − 1) (36hxn + 9h2 s + 26h2 + 10xn2 + 4hsxn) fn+1 − 1 480h5(s − 2) (32hxn + 8h2 s + 19h2 + 10xn2 + 4hsxn) fn+2 + 7h2 + 8hxn + 2xn2 24h5 s ∗ (s4 − 10s3 + 35s2 − 50s + 24 fn+s + 1 720h5(s − 3) (28hxn + 7h2 s + 14h2 + 10xn2 + 4hsxn) fn+3 − 1 2880h5(s − 4) (24hxn + 6h2 s + 11h2 + 10xn2 + 4hsxn) fn+4 a7 = 10h + 5xn + hs 5040h5 s fn − 9h + 5xn + hs 1260h5(s − 1) fn+1 + 8h + 5xn + hs 840h5(s − 2) fn+2 − 2h + xn 42h5 s(s4 − 10s3 + 35s2 − 50s + 24) fn+s − 7h + 5xn + hs 12605(s − 3) fn+3 + 6h + 5xn + hs 5040h5(s − 4) fn+4 a8 = − 1 8064h5 s(s4 − 10s3 + 35s2 − 50s + 24) (24 fn − 24 fn+s − 50s fn + 96s fn+1 75 J. O. Kuboye et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 69–76 76 −72s fn+2 + 32s fn+3 − 6s fn+4 + 35s 2 fn − 10s 3 fn − 104s 2 fn+1 + s 4 fn + 36s 3 fn+1 +114s2 fn+2 +11s2 fn+4 − 4s 4 fn+3 − 6s 3 fn+4 + s 4 fn+4 ) Appendix B E =  20480 371589120 90112 371589120 28672 371589120 315392 371589120 430080 371589120 585111 371589120 665797 371589120 236334 371589120 4096 10321920 24576 10321920 14336 10321920 86016 10321920 129024 10321920 262947 10321920 346473 10321920 132678 10321920 −20480 6881280 −155648 6881280 −229376 6881280 602112 6881280 1290240 6881280 496249 6881280 2271083 6881280 1180914 6881280 4096 1270080 32768 1270080 57344 1270080 114688 1270080 344064 1270080 487059 1270080 351161 1270080 64806 1270080 −20480 30965760 −188416 30965760 −501760 30965760 143360 30965760 2795520 30965760 4595169 30965760 2242979 30965760 490434 30965760 20480 289013760 221184 289013760 888832 289013760 1634304 289013760 1290240 289013760 −837513 3870720 −693915 289013760 −200466 289013760  F =  −163840 371589120 −630784 371589120 172032 371589120 1576960 371589120 1720320 371589120 −1170222 371589120 −665797 371589120 32768 10321920 172032 10321920 86016 10321920 −430080 10321920 −516096 10321920 525894 10321920 346473 10321920 −163840 6881280 −1089536 6881280 −1376256 6881280 3010560 6881280 5160960 6881280 992498 6881280 2271083 6881280 32768 1270080 229376 1270080 344064 1270080 −573440 1270080 −3176256 1270080 974118 1270080 351161 1270080 −163840 30965760 −1318912 30965760 −3010560 30965760 716800 30965760 11182080 30965760 9190338 30965760 2242979 30965760 163840 289013760 1548288 289013760 752645332992 289013760 8171520 289013760 5160960 289013760 −1675026 289013760 −693915 289013760  76