J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 Journal of the Nigerian Society of Physical Sciences Numerical Solution of Stiff and Oscillatory Problems using Third Derivative Trigonometrically Fitted Block Method M. Kidaa, S. Adamub,∗, O. O. Adurojac, T. P. Pantuvod aDepartment of Mathematics, Modibbo Adama University Yola, Adamawa, Nigeria bDepartment of Mathematics, Nigerian Army University Biu, Borno, Nigeria cDepartment of Mathematics, Osun State College of Education Ilesha, Osun, Nigeria dMathematics and Statistics Department, Federal University Wukari, Taraba, Nigeria Abstract This paper considered the formulation of continuous third derivative trigonometrically fitted block method for the solution of stiff and oscillatory problems. The development of the technique involved the interpolation and collocation of the approximate solution which is the combination of polynomial and trigonometric functions. Solving for the unknown parameters and substituting the results into the approximate solution yielded a continuous linear multistep method, which is evaluated at some selected grid points where two cases were considered at equal intervals to give the discrete schemes which are implemented in block form. The blocks are convergent and stable. Numerical experiments show that the methods compete favorably with existing method and efficient for the solution of stiff and oscillatory problems. DOI:10.46481/jnsps.2022.271 Keywords: Third derivative, Trigonometrically fitted block method, Oscillatory problems, Stability, Convergence Article History : Received: 22 June 2021 Received in revised form: 27 January 2022 Accepted for publication: 28 January 2022 Published: 28 February 2022 c©2022 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction In science and engineering, mathematical models are devel- oped to help in the studying of physical phenomena, these mod- els often yield equations that contain some derivatives of an un- known function of one or several variables, such equations are called Differential Equations (DEs) [1, 2]. Interestingly, sys- tems described by differential equations are so complex, or the systems that they describe are so large, that a purely analyti- cal solution to the equations is not tractable, hence the need for numerical approximation [3]. ∗Corresponding author tel. no: +2348069405058 Email address: malgwisa@gmail.com (S. Adamu ) Block method is formulated in terms of LMM, It preserves the traditional advantage of one step methods, of being self- starting and permitting easy change of step length [1, 4]. The advantages of block methods over predictor-corrector methods lies in the fact that they are less expensive in terms of number of functions evaluation, it is capable of giving evaluations at dif- ferent grid points, without overlapping as done in the predictor- corrector method and it generates simultaneous solutions at all grid points [1, 4, 5]. Despite the success recorded by linear multistep method (LMM) in the numerical solution of initial value problems, most of the approaches failed when the problems are stiff. Adoption of trigonometrically fitted approximate solution as basis func- tions have been effective in handling this setback, but their de- 34 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 35 velopments are tedious. In this paper, we consider numerical procedures for approximating the solution of stiff and oscilla- tory problems of the form y′ = f (x, y, ), y (x0) = η0. (1) These type of differential equations are known to be highly os- cillatory and some problems have special properties such as dis- continuity and stiffness, hence it is quite difficult to get their numerical solutions accurately [6, 2]. They are numerically un- stable, unless the step size is taken to be extremely small [7]. A stiff system is one involving rapidly changing components to- gether with slowly changing ones, stiff systems are sometimes referred to as systems with large Lipstchitz constant [7, 8]. The system of (1) is stiff if the eigenvalues of the matrix have neg- ative real parts at every time x and varies greatly in magnitude [9, 10]. Oscillation is the repetitive variation, typically in time, of some measure about a central value (often a point of equilib- rium) or between two or more different states [7]. This type of problem arises in different fields of science and engineering, which includes quantum mechanics, celestial me- chanics, molecular dynamics, quantum chemistry, astrophysics, electronics and semi-discretizations of wave equation [5, 11]. Numerous numerical methods have been derived for ap- proximating the solutions of (1), some of which are Linear Mul- tistep Methods (LMMs) [12, 13, 4, 14], Block method [5, 15, 1, 16, 17, 18, 19, 20], Trigonometrically fitted methods [21, 22, 23, 6, 24]. Other methods for the solution of Stiff and oscil- latory problems includes [2, 9, 11]. Modeled problems in en- gineering and sciences leads to complex problems such as stiff and oscillatory problems. It was observed that LMM fails to converge at the state of oscillatory as a result has poor stability properties [7]. Second derivative formulas with trigonometric basis functions seems to solve such problems very well and has a better and stability properties [1, 17]. [25, 26] develop an algorithm for the solution of Stiff Initial Value Problems. Furthermore, [27, 28] considers an implicit r−point block backward differentiation formula (BBDF) and block method that generates two values simultaneously respectively for solv- ing stiff IVPs of ODEs. The r−point block method simultane- ously produces r−new values at the time discretization points. The total number of steps to complete the integration by the methods are reduced but not sufficiently and the computation time for the method can still be reduced. [29] developed a block integrator for the solution of stiff and oscillatory first order IVPs of ODEs by means of collocation and interpolation of the com- bination of power series and exponential function to generate a continuous implicit linear multistep method. Despite the success achieved by the above methods, it in- creases the dimension of the resulting systems of first order by the order of the differential equation, hence it wastes both the computer and human effort. Polynomial approximate solution has been reported not to be efficient in handling stiff oscillatory problems, trigonometrically fitted methods handle oscillatory problems effectively, but the application of higher derivatives methods for the solution of first order oscillatory initial value problems has not given much attention. In this paper Continuous Trigonometrically Fitted Third Deriva- tive Method (CTFTDM) is constructed which provides a dis- crete method for direct solution of first order initial value prob- lems. The coefficients of the TFTDM are functions of the fre- quency and the step-size, hence the solutions that will be pro- vided by the methods will be highly accurate if (1) has periodic solutions with the unknown frequencies. This new method con- siders the application of higher derivative and does not waste both the computer and human effort which is efficient in han- dling stiff oscillatory problems. Trigonometrically fitted meth- ods are powerful tool in handling stiff/oscillatory problems, or- der 2 methods developed by combining polynomial and trigono- metric approximate solution and considering higher derivatives, conveniently solved stiff/oscillatory problems and it perform better than other methods of the same order. 1.1. Preliminaries Definition 1.1. Stiff [8]: The system ẏ = Ay + φ(x) where A is an m × m matrix and φ(x) an m dimensional vector, is said to be stiff if Reλi < 0, i = 1, 2, . . . , m and maxi=1...m |Reλi| mini=1,...,m |Reλi| ≥ 0 where λi, i = 1, 2, . . . , m are the eigenvalues of A. The ratio S = maxi=1...m |Reλi| mini=1,...,m |Reλi| ≥ 0 is called the stiffness ratio. Theorem 1.2. Consistency [30]: A block method is said to be consistent if it has order p ≥ 1. Theorem 1.3. Zero-stable [30]: A block method is said to be zero stable if as h → 0, the roots r j, j = 1(1)k of the first char- acteristic polynomials ρ(r) = 0 that is ρ(r) = det [∑ A(0)Rk−1 ] = 0 satisfying |R| ≤ 1, must be simple. Theorem 1.4. Convergent [31]: Consistent and zero stability are sufficient conditions for a block method to be convergent. 2. Methodology We considered the approximate solution y (x) = 2∑ n=0 αn x n + 2∑ n=1 βn sin nωx + 2∑ n=1 cn cos nωx (2) 35 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 36 where ω = 2nπ T , T is the period of the oscillation, α′s, β′s and c′s are parameters to be determined. Interpolating (2) at point xn, collocating the first derivative of (2) at point xn, xn+u, xn+v, collocating the second derivative of (2) at point xn+u, xn+v, and collocating the third derivative of (2) at point xn+v for u and v, u < v to obtain the system XA = U (3) where A = [ α1 β1 β2 β3 c1 c2 d2 ]T U = [ yn fn fn+u fn+v gn+u gn+v ln+v ]T X =  1 xn x2n sin ωxn sin 2ωxn cos ωxn cos 2ωxn 0 1 2xn ω cos ωxn+v 2ω cos 2ωxn −ω sin ωxn −2ω sin 2ωxn 0 1 2xn+u ω cos ωxn+u 2ω cos 2ωxn+u −ω sin ωxn+u −2ω sin 2ωxn+u 0 1 2xn+v ω cos ωωxn+v 2ω cos 2ωxn+v −ω sin ωxn+v −2ω sin 2ωxn+v 0 0 2 −ω2 sin ωxn+u −4ω2 sin 2ωxn+u −ω2 cos ωxn+u −4ω2 cos 2ωxn+u 0 0 2 −ω2 sin ωxn+v −4ω2 sin 2ωxn+v −ω2 cos ωxn+v −4ω2 cos 2ωxn+v 0 0 0 −ω3 cos ωxn+v −8ω3 cos 2ωxn+v ω3 sin ωxn+v 8ω3 sin 2ωxn+v  (3) gives system of nonlinear equations which is then solved using Newton Raphson’s method. Now considering two cases of third derivative trigonometrically fitted methods at equal intervals as specified below Case I: u = 12 , v = 1 Case II: u = 1, v = 2 2.1. Block Method for Case I Considering u = 12 , v = 1 yn+t = α0 (t) yn + β0 (t) fn + β1 (t) fn+ 12 + β2 (t) fn+1 + γn+1 (t) gn+ 12 (4) +γn+2 (t) gn+2 + [ ln+2 (t) (ζn+2) ] Evaluating (4) at t = 12 yields yn+ 12 = α01yn + β01 fn + β11 fn+ 12 + β21 fn+1 + γ11gn+ 12 (5) +γ21gn+1 + ζ11ln+1 where α01 = 1 β01 = 1 4  8 cos hω + 60 cos 2hω + 104 cos 12 hω− 94 cos 3 2 hω− 10 cos 5 2 hω + 30h 2ω2 +6h2ω2 cos hω− 33h2ω2 cos 12 hω− 3h 2ω2 cos 32 hω + 8hω sin hω +8hω sin 2hω− 56hω sin 12 hω + 6hω sin 3 2 hω− 2hω sin 5 2 hω− 68  ω [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] β11 = 1 4  164 cos hω− 44 cos 2hω− 4 cos 3hω + 20 cos 12 hω− 42 cos 3 2 hω +22 cos 52 hω + 18h 2ω2 + 6h2ω2 cos 2hω− 30h2ω2 cos 12 hω +9h2ω2 cos 32 hω− 3h 2ω2 cos 52 hω + 32hω sin hω− 40hω sin 2hω +28hω sin 12 hω + 6hω sin 3 2 hω + 10hω sin 5 2 hω− 116  ω [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] β21 = 1 4  −172 cos hω− 16 cos 2hω + 4 cos 3hω− 124 cos 12 hω + 136 cos 3 2 hω −12 cos 52 hω + 4h 2ω2 + 2h2ω2 cos hω− 2h2ω2 cos 2hω− h2ω2 cos 12 hω −4h2ω2 cos 32 hω + h 2ω2 cos 52 hω− 16hω sin hω− 12hω sin 2hω −56hω sin 12 hω + 42hω sin 3 2 hω + 2hω sin 5 2 hω + 184  ω [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] 36 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 37 γ11 = 1 4  54 sin hω + 28 sin 2hω + 6 sin 3hω + 12 sin 12 hω− 66 sin 3 2 hω− 14 sin 5 2 hω −138hω + 6h2ω2 sin 2hω + 18h2ω2 sin 12 hω− 9h 2ω2 sin 32 hω −3h2ω2 sin 52 hω− 6hω cos hω + 18hω cos 2hω− 2hω cos 3hω +156hω cos 12 hω− 22hω cos 3 2 hω− 6hω cos 5 2 hω  ω2 [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] γ21 = 1 8  −108 sin hω− 56 sin 2hω− 12 sin 3hω− 24 sin 12 hω + 132 sin 3 2 hω +24hω + 12h2ω2 sin hω + 30h2ω2 sin 2hω + 6h2ω2 sin 12 hω− 45h 2ω2 sin 32 hω −3h2ω2 sin 52 hω + 192hω cos hω + 40hω cos 2hω− 84hω cos 1 2 hω +28 sin 52 hω− 186hω cos 3 2 hω + 14hω cos 5 2 hω  ω2 [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] ζ11 = 1 8  −156 cos hω− 32 cos 2hω− 4 cos 3hω− 144 cos 12 hω + 136 cos 3 2 hω +2h2ω2 + 16h2ω2 cos hω + 14h2ω2 cos 2hω− 2h2ω2 cos 12 hω −h2ω2 cos 52 hω− 144hω sin hω− 24hω sin 2hω + 36hω sin 1 2 hω +8 cos 52 hω− 29h 2ω2 cos 32 hω + 126hω sin 3 2 hω− 6hω sin 5 2 hω + 192  ω3 [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] Evaluating (4) at t = 1 yn+1 = α02yn + β02 fn + β12 fn+ 12 + β22 fn+1 + γ12gn+ 12 (6) +γ22gn+1 + ζ12ln+1 Where α02 = 1 β02 = − 1 2  32 cos hω− 30 cos 2hω− 34 cos 12 hω + 29 cos 3 2 hω + 5 cos 5 2 hω− 20h 2ω2 +18h2ω2 cos 12 hω + 2h 2ω2 cos 32 hω− 16hω sin hω− 4hω sin 2hω +70hω sin 12 hω− 9hω sin 3 2 hω + hω sin 5 2 hω− 2  ω  12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω −hω cos 52 hω  β12 =  65 cos hω− 20 cos 2hω− cos 3hω + 2 cos 12 hω− 9 cos 3 2 hω +7 cos 52 hω + 6h 2ω2 + 4h2ω2 cos hω + 2h2ω2 cos 2hω +3h2ω2 cos 32 hω− h 2ω2 cos 52 hω + 8hω sin hω− 16hω sin 2hω −14h2ω2 cos 12 hω + 28hω sin 1 2 hω + 4hω sin 5 2 hω− 44  ω  12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω −hω cos 52 hω  β22 = − 1 2  98 cos hω− 10 cos 2hω− 2 cos 3hω + 38 cos 12 hω− 47 cos 3 2 hω +9 cos 52 hω− 20h 2ω2 + 18h2ω2 cos 12 hω + 2h 2ω2 cos 32 hω +8hω sin hω + 16hω sin 2hω + 70hω sin 12 hω− 45hω sin 3 2 hω −3hω sin 52 hω− 86  ω  12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω +26hω + 4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω +hω cos 32 hω− hω cos 5 2 hω  γ12 = 1 2  39 sin hω− 24 sin 2hω + 3 sin 3hω− 12 sin 12 hω− 6 sin 3 2 hω +6 sin 52 hω− 26hω + 8h 2ω2 sin hω + 2h2ω2 sin 2hω + 4h2ω2 sin 12 hω −6h2ω2 sin 32 hω− 2h 2ω2 sin 52 hω + 17hω cos hω + 10hω cos 2hω −hω cos 3hω + 8hω cos 12 hω− 8hω cos 5 2 hω  ω2  12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω +26hω− hω cos 52 hω  37 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 38 γ22 = 1 2  −39 sin hω + 24 sin 2hω− 3 sin 3hω + 12 sin 12 hω + 6 sin 3 2 hω −6 sin 52 hω− 16hω + 4h 2ω2 sin hω + 8h2ω2 sin 2hω −4h2ω2 sin 12 hω− 12h 2ω2 sin 32 hω + 16hω cos hω + 28hω cos 1 2 hω −38hω cos 32 hω + 10hω cos 5 2 hω  ω2  12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω +26hω + 4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω −hω cos 52 hω  ζ12 = 1 2  −63 cos hω + 10 cos 2hω− cos 3hω− 24 cos 12 hω + 28 cos 3 2 hω −4 cos 52 hω + 8h 2ω2 + 4h2ω2 cos hω + 4h2ω2 cos 2hω −8h2ω2 cos 12 hω− 8h 2ω2 cos 32 hω− 24hω sin hω− 4hω sin 2hω −24hω sin 12 hω + 36hω sin 3 2 hω− 4hω sin 5 2 hω + 54  ω3  12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω +26hω + 4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω +hω cos 32 hω− hω cos 5 2 hω  writing (5) and (6) in discrete block method[ 1 0 0 1 ] [ yn+ 12 yn+1 ] = [ 0 1 0 1 ] [ yn−1 yn ] + [ 0 β01 0 β02 ] [ fn−1 fn ] (7) + [ β11 β21 β12 β22 ] [ fn+ 12 fn+1 ] + [ γ11 γ12 γ21 γ22 ] [ gn+ 12 gn+1 ] + [ 0 ζ11 0 ζ21 ] [ ln−1 ln+1 ] 2.2. Analysis of the Stability Properties for Case I 2.2.1. Order and error constant Evaluating each row of (5) and (6) in Taylor series about xn gives L [ y (x) ; h ] = yn+ 12 − yn − hy ′ n −β01 fn −β11 fn+ 12 −β21 fn+1 −γ11gn+ 12 −γ21gn+1 −γ11ln+1 = 0 therefore, the block methods are of order p = (2, 2)T with the following error constants k = − 1 96  −1872 cos hω− 384 cos 2hω− 48 cos 3hω− 1728 cos 12 hω + 1632 cos 3 2 hω +96 cos 52 hω + 516h 2ω2 + 50h4ω4 + 852h2ω2 cos hω + 540h2ω2 cos 2hω +12h2ω2 cos 3hω + 16h4ω4 cos hω− 10h4ω4 cos 2hω + 24h3ω3 sin hω +212h3ω3 sin 2hω− 588h2ω2 cos 12 hω− 1338h 2ω2 cos 32 hω + 6h 2ω2 cos 52 hω −38h4ω4 cos 12 hω− 23h 4ω4 cos 32 hω + 5h 4ω4 cos 52 hω− 216h 3ω3 sin 12 hω −180h3ω3 sin 32 hω− 28h 3ω3 sin 52 hω− 2376hω sin hω− 624hω sin 2hω −72hω sin 3hω + 288hω sin 12 hω + 2304hω sin 3 2 hω + 96hω sin 5 2 hω + 2304  ω3 [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] k = 1 24  756 cos hω− 120 cos 2hω + 12 cos 3hω + 288 cos 12 hω− 336 cos 3 2 hω + 48 cos 5 2 hω −132h2ω2 − 34h4ω4 + 51h2ω2 cos hω− 108h2ω2 cos 2hω− 3h2ω2 cos 3hω +4h4ω4 cos hω + 2h4ω4 cos 2hω− 24h3ω3 sin hω− 52h3ω3 sin 2hω −66h2ω2 cos 12 hω + 297h 2ω2 cos 32 hω− 39h 2ω2 cos 52 hω + 22h 4ω4 cos 12 hω +7h4ω4 cos 32 hω− h 4ω4 cos 52 hω + 192h 3ω3 sin 12 hω + 18h 3ω3 sin 32 hω +2h3ω3 sin 52 hω + 522hω sin hω− 96hω sin 2hω + 18hω sin 3hω + 216hω sin 1 2 hω −468hω sin 32 hω + 84hω sin 5 2 hω− 648  ω3 [ 12 sin hω− 22 sin 2hω− 42 sin 12 hω + 27 sin 3 2 hω + 5 sin 5 2 hω + 26hω +4hω cos hω + 2hω cos 2hω− 32hω cos 12 hω + hω cos 3 2 hω− hω cos 5 2 hω ] 2.2.2. Zero stability of the block − ρ(λ) = det [ λA(1) − A(0) ] = 0 (8) 38 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 39 where A(0), A(1) are from (7) λ [ 1 0 0 1 ] − [ 0 1 0 1 ] = det [ λ −1 0 λ− 1 ] ,λ2 −λ = λ (λ− 1) = 0 Thus, solving for λ λ (λ− 1) = 0 which implies that λ1 = 0,λ2 = 1. Hence, using theorem 1, theorem 2, theorem 3, the block method (7) is zero stable and also consistent as its order p = [2, 2]T > 1, thus it is convergent. 2.2.3. Linear stability Applying the test equation y(k) = λ(k)yn to yield yn+1 = µ (z) ym, µ (z) is the amplification equation given by µ (z) = − ( A(1) − zβ(1) − zγ(1) − zζ(1) )−1 ( A(0) + z2β(0) + z3γ(0) ) (9) where A(0) = [ 0 1 0 1 ] , β(0) = [ 0 β01 0 β02 ] , γ(0) = [ 0 ζ11 0 ζ21 ] , A(1) = [ 1 0 0 1 ] , β(1) = [ β11 β21 β12 β22 ] and γ(1) = [ γ11 γ12 γ21 γ22 ] the matrix µ (z) has eigenvalues (0, 0, ...,ξk) where ξk is called the stability function which is the rational function with coefficients. µ (z) =  72z − 234z2h + 9zh3 − 102z2h2 + 27z2h4 + 126zh −60h3ω4 − 64h4ω4 + 210zh2ω2 + 210zh3ω2 − 108zh4ω4 −136zh5ω4 + 180z2h3ω2 + 150z2h4ω2 + 36z2h5ω4 − 24z2h6ω4 [ 72z − 180z2h + 42z3h2 + 45z4h3 − 60h3ω4 + 210zh2ω2 −112zh4ω4 + 180z2h3ω2 − 135z3h4ω2 + 42z2h5ω4 ] , 0 2.2.4. Region of absolute stability The region of absolute stability of the block method (7) is shown in Figure 1. Figure 1. Region of Absolute Stability for Case I 2.3. Block Method for Case II We considered u = 1, v = 2 yn+t = α0 (t) yn + β0 (t) fn + β1 (t) fn+1 + β2 (t) fn+2 + γn+1 (t) gn+1 (10) +γn+2 (t) gn+2 + ζn+2 (t) ln+2 Evaluating (4) at t = 1 yields yn+1 = α01yn + [ β01 fn + β11 fn+1 + β21 fn+2 ] + [ γ11gn+1 + γ21gn+2 ] + ζ11ln+2 (11) 39 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 40 where α01 = 1 β01 = 1 2  −52 cos hω− 4 cos 2hω + 47 cos 3hω− 30 cos 4hω +5 cos 5hω− 60h2ω2 + 66h2ω2 cos hω− 12h2ω2 cos 2hω +6h2ω2 cos 3hω + 56hω sin hω− 8hω sin 2hω− 6hω sin 3hω −8hω sin 4hω + 2hω sin 5hω + 34  ω  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  β11 = 1 2  −10 cos hω− 82 cos 2hω + 21 cos 3hω + 22 cos 4hω −11 cos 5hω + 2 cos 6hω− 36h2ω2 + 60h2ω2 cos hω −18h2ω2 cos 3hω− 12h2ω2 cos 4hω + 6h2ω2 cos 5hω −28hω sin hω− 32hω sin 2hω− 6hω sin 3hω + 40hω sin 4hω −10hω sin 5hω + 58  ω  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  β21 =  31 cos hω + 43 cos 2hω− 34 cos 3hω + 4 cos 4hω + 3 cos 5hω −cos 6hω− 4h2ω2 + h2ω2 cos hω− 2h2ω2 cos 2hω +4h2ω2 cos 3hω + 2h2ω2 cos 4hω− h2ω2 cos 5hω +28hω sin hω + 8hω sin 2hω− 21hω sin 3hω +6hω sin 4hω− hω sin 5hω− 46  ω  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  γ11 = − 1 2  6 sin hω + 27 sin 2hω− 33 sin 3hω + 14 sin 4hω− 7 sin 5hω +3 sin 6hω− 138hω + 36h2ω2 sin hω− 18h2ω2 sin 3hω +12h2ω2 sin 4hω− 6h2ω2 sin 5hω + 156hω cos hω −6hω cos 2hω− 22hω cos 3hω + 18hω cos 4hω −6hω cos 5hω− 2hω cos 6hω  ω2  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  γ21 = − 1 2  −6 sin hω− 27 sin 2hω + 33 sin 3hω− 14 sin 4hω + 7 sin 5hω −3 sin 6hω + 12hω + 6h2ω2 sin hω + 12h2ω2 sin 2hω −45h2ω2 sin 3hω + 30h2ω2 sin 4hω− 3h2ω2 sin 5hω− 42hω cos hω +96hω cos 2hω− 93hω cos 3hω + 20hω cos 4hω + 7hω cos 5hω  ω2  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  ζ11 = 1 2  36 cos hω + 39 cos 2hω− 34 cos 3hω + 8 cos 4hω− 2 cos 5hω + cos 6hω− 2h2ω2 + 2h2ω2 cos hω− 16h2ω2 cos 2hω +29h2ω2 cos 3hω− 14h2ω2 cos 4hω + h2ω2 cos 5hω −18hω sin hω + 72hω sin 2hω− 63hω sin 3hω + 12hω sin 4hω +3hω sin 5hω− 48  ω3  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  Evaluating (4) at t = 2 yields yn+2 = α02yn + [ β02 fn + β12 fn+1 + β22 fn+2 ] + [ γ12gn+1 + γ22gn+2 ] + ζ12ln+2 (12) 40 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 41 where α01 = 1 β02 = 1 2  −34 cos hω + 32 cos 2hω + 29 cos 3hω− 30 cos 4hω +5 cos 5hω− 80h2ω2 + 72h2ω2 cos hω + 8h2ω2 cos 3hω +140hω sin hω− 32hω sin 2hω− 18hω sin 3hω −8hω sin 4hω + 2hω sin 5hω− 2  ω  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  β12 =  −2 cos hω− 65 cos 2hω + 9 cos 3hω + 20 cos 4hω −7 cos 5hω + cos 6hω− 24h2ω2 + 56h2ω2 cos hω −16h2ω2 cos 2hω− 12h2ω2 cos 3hω− 8h2ω2 cos 4hω +4h2ω2 cos 5hω− 56hω sin hω− 16hω sin 2hω +32hω sin 4hω− 8hω sin 5hω + 44  ω  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  β22 = 1 2  38 cos hω + 98 cos 2hω− 47 cos 3hω− 10 cos 4hω +9 cos 5hω− 2 cos 6hω− 80h2ω2 + 72h2ω2 cos hω +8h2ω2 cos 3hω + 140hω sin hω + 16hω sin 2hω −90hω sin 3hω + 32hω sin 4hω− 6hω sin 5hω− 86  ω  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  γ12 = − 1 2  −12 sin hω + 39 sin 2hω− 6 sin 3hω− 24 sin 4hω + 6 sin 5hω +3 sin 6hω− 52hω + 16h2ω2 sin hω + 32h2ω2 sin 2hω −24h2ω2 sin 3hω + 8h2ω2 sin 4hω− 8h2ω2 sin 5hω +16hω cos hω + 34hω cos 2hω + 20hω cos 4hω −16hω cos 5hω− 2hω cos 6hω  ω2  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  γ22 = 1 2  −12 sin hω + 39 sin 2hω− 6 sin 3hω− 24 sin 4hω +6 sin 5hω + 3 sin 6hω + 32hω + 16h2ω2 sin hω −16h2ω2 sin 2hω + 48h2ω2 sin 3hω− 32h2ω2 sin 4hω −56hω cos hω− 32hω cos 2hω + 76hω cos 3hω− 20hω cos 5hω  ω2  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  ζ12 = 1 2  24 cos hω + 63 cos 2hω− 28 cos 3hω− 10 cos 4hω +4 cos 5hω + cos 6hω− 32h2ω2 + 32h2ω2 cos hω −16h2ω2 cos 2hω + 32h2ω2 cos 3hω− 16h2ω2 cos 4hω +48hω sin hω + 48hω sin 2hω− 72hω sin 3hω + 8hω sin 4hω +8hω sin 5hω− 54  ω3  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  writing (11) and (12) in discrete block method[ 1 0 0 1 ] [ yn+1 yn+2 ] = [ 0 1 0 1 ] [ yn−1 yn ] + [ 0 β01 0 β02 ] [ fn−1 fn ] (13) + [ β11 β21 β12 β22 ] [ fn+1 fn+2 ] + [ γ11 γ12 γ21 γ22 ] [ gn+1 gn+2 ] + [ 0 ζ11 0 ζ21 ] [ ln−1 ln+1 ] 41 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 42 2.4. Analysis of the Stability Properties for Case II 2.4.1. Order and error constant Evaluating each row of (11) and (12) in Taylor series about xn gives L [ y (x) ; h ] = yn+1 − yn −β01 fn −β11 fn+1 −β21 fn+2 −γ11gn+1 −γ21gn+2 − ζ11ln+2 = 0 therefore, the methods are of order p = (2, 2)T with the following error constants k = − 1 12  216 cos hω + 234 cos 2hω− 204 cos 3hω + 48 cos 4hω −12 cos 5hω + 6 cos 6hω− 258h2ω2 − 100h4ω4 + 294h2ω2 cos hω −426h2ω2 cos 2hω + 669h2ω2 cos 3hω− 270h2ω2 cos 4hω +76h4ω4 cos hω− 3h2ω2 cos 5hω− 6h2ω2 cos 6hω− 32h4ω4 cos 2hω +46h4ω4 cos 3hω + 20h4ω4 cos 4hω− 10h4ω4 cos 5hω + 216h3ω3 sin hω −24h3ω3 sin 2hω + 180h3ω3 sin 3hω− 212h3ω3 sin 4hω + 28h3ω3 sin 5hω −72hω sin hω + 594hω sin 2hω− 576hω sin 3hω + 156hω sin 4hω −24hω sin 5hω + 18hω sin 6hω− 288  ω3  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω− 2hω cos 3hω −4hω cos 4hω + 2hω cos 5hω  k = − 1 6  72 cos hω + 189 cos 2hω− 84 cos 3hω− 30 cos 4hω + 12 cos 5hω +3 cos 6hω− 132h2ω2 − 136h4ω4 − 66h2ω2 cos hω + 51h2ω2 cos 2hω +297h2ω2 cos 3hω− 108h2ω2 cos 4hω + 88h4ω4 cos hω− 39h2ω2 cos 5hω −3h2ω2 cos 6hω + 16h4ω4 cos 2hω + 28h4ω4 cos 3hω + 8h4ω4 cos 4hω −4h4ω4 cos 5hω + 384h3ω3 sin hω− 48h3ω3 sin 2hω + 36h3ω3 sin 3hω −104h3ω3 sin 4hω + 4h3ω3 sin 5hω + 108hω sin hω + 261hω sin 2hω −234hω sin 3hω− 48hω sin 4hω + 42hω sin 5hω + 9hω sin 6hω− 162  ω3  42 sin hω− 12 sin 2hω− 27 sin 3hω + 22 sin 4hω −5 sin 5hω− 52hω + 64hω cos hω− 8hω cos 2hω −2hω cos 3hω− 4hω cos 4hω + 2hω cos 5hω  2.4.2. Zero stability of the block − ρ(λ) = det [ λA(1) − A(0) ] = 0 (14) where A(0), A(1) are from (13) λ [ 1 0 0 1 ] − [ 0 1 0 1 ] = det [ λ −1 0 λ− 1 ] ,λ2 −λ = λ (λ− 1) = 0 Thus, solving for λ λ (λ− 1) = 0 which implies that λ1 = 0,λ2 = 1 Hence, using theorem 1, theorem 2, theorem 3, the block method (13) is zero stable and also consistent as its order p = [2, 2]T > 1, thus it is convergent. 2.4.3. Linear stability Applying the test equation y(k) = λ(k)yn to yield yn+1 = µ (z) ym, µ (z) is the amplification equation given by µ (z) = − ( A(1) − zβ(1) − z2γ(1) − z3ζ(1) )−1 ( A(0) + zβ(0) + z2γ(0) ) (15) where A(0) = [ 0 1 0 1 ] , β(0) = [ 0 β01 0 β02 ] , γ(0) = [ 0 ζ11 0 ζ21 ] , A(1) = [ 1 0 0 1 ] , β(1) = [ β11 β21 β12 β22 ] and γ(1) = [ γ11 γ12 γ21 γ22 ] the matrix µ (z) has eigenvalues (0, 0, ...,ξk) where ξk is called the stability function which is the rational function with coefficients. µ (z) = [ 1364h5z4ω4 − 2816h4z2ω4 − 180h3z5 + 280h3z4ω2 +960h3ω4 + 270h2z3 − 420h2z2ω2 − 246hz4 + 369z2 ]  −288h5z4ω4 + 288h5z3ω4 + 560h4z2ω4 − 1260h3z4ω2 +360h3z3ω2 − 960h3ω4 + 150h2z5 + 2100h2z2ω2 +162hz4 − 342hz3 − 225z2  , 0 42 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 43 Table 1. Result for Example I h 3BEBDF MAXE 3BBDF MAXE Error Case I Time Error Case II Time 10−2 1.68449e(−001) 6.62694e(−002) 3.6891e(−002) 0.178s 2.9440e(−002) 0.178s 10−3 5.14997e(−002) 7.44768e(−002) 3.6837e(−003) 0.178s 8.3499e(−003) 0.178s 10−4 6.95725e(−003) 8.45376e(−003) 2.6596e(−004) 0.179s 3.8627e(−004) 0.179s 10−5 7.16727e(−004) 8.53717e(−004) 2.6585e(−005) 0.179s 1.0017e(−005) 0.179s 2.4.4. Region of absolute stability The region of absolute stability of the block method (13) is shown in Figure 2 Figure 2. Region of Absolute Stability for Case II 3. Numerical Experiments Evaluating the performance of the new block method on some challenging stiff and oscillatory problems which appeared in literature and compared the results with solution from some methods. The following notations were used in table: h Step size MAXE Maximum Error 3BBDF 3-point block backward differentiation formula 3BE BDF 3-point block extended backward differentiation formula Example I: Consider a linear oscillatory problem in the interval 0 ≤ x ≤ 10( y′1 = 9y1 + 24y2 + 5 cos x − 1 3 sin x y′2 = −24y1 − 5y2 − 9 cos x + 1 3 sin x ) , ( y1 (0) y2 (0) ) = ( 4 3 2 3 ) The exact solution is given as ( y1 (x) y2 (x) ) = ( 4e−x − 3e−1000x −2e−x + 3e−1000x ) Source: [25] Table 1 shows that the new methods performed accurately and approximate better than the result of [25]. For instance, at h = 10−5, the error in [25] are 7.16727e(−004) and 8.53717e(−004) respectively, while the error in case 1 and case 2 are 2.6585e(−05) and 1.0017e(−05) respectively. The comparison of error in Figure 3 clearly shows that, the new methods converge better than the existing method. 43 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 44 Figure 3. Error for Example I Table 2. Results for Example II h 3BEBDF MAXE 3BBDF MAXE CASE I Time CASE II Time 10−2 8.35359e(−002) 1.08041e(−002) 1.7729e(−05) 0.3511s 1.7729e(−05) 0.3511s 10−3 9.10218e(−003) 1.08962e(−002) 1.0051e(−04) 0.3511s 2.1635e(−06) 0.3511s 10−4 9.18073e(−004) 1.09230e(−003) 5.0175e(−05) 0.3511s 2.0646e(−07) 0.3511s 10−5 9.18862e(−005) 1.09257e(−004) 4.1886e(−06) 0.3510s 7.7533e(−08) 0.3510s 10−6 9.18939e(−006) 1.09259e(−005) 4.3258e(−07) 0.3510s 7.4327e(−09) 0.3510s Example II: Consider the linear stiff system in the interval 0 ≤ x ≤ 10( y′1 y′2 ) = ( −100y1 + 9.901y2 0.1y1 − y2 ) , ( y1 (0) y2 (0) ) = ( 1 10 ) The eigenvalues of the Jacobian matrix of the system are λ1= −0.99 and λ2 = −100.01. The exact solution is given as( y1 (x) y2 (x) ) = ( exp (−0.99x) 10 exp (−0.99x) ) Source: [25] Table 2 shows that the new methods performed accurately and approximate better than the result of [25]. At h = 10−6, the error of [25] are 9.18939e(−006) and 1.09259e(−005) , while the error in case 1 and case 2 are 4.3258e(−07) and 7.4327e(−09) respectively. This clearly shows that the new methods performed better than of [25].The comparison of error in Figure 4 clearly shows that, the new methods converge better than the existing method. Example III: Consider a linear oscillatory problem in the interval 0 ≤ x ≤ 10( y′1 y′2 ) = ( −2y1 + y2 + 2 sin x 998y1 − 999y2 + 999 (cos x − sin x) ) , ( y1 (0) y2 (0) ) = ( 2 3 ) The exact solution is given as ( y1 (x) y2 (x) ) = ( 2 exp (−x) + sin x 2 exp (−x) + cos x ) Source: [26] Table 3 shows that the new methods performed accurately and approximate better than the result of [26]. From the result obtained for example III as shown in the Table 3, it is evident that the class of third derivative trigonometrically fitted method 44 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 45 Figure 4. Error for Example II Table 3. Result for Example III h MAXE in [26] CASE I Time CASE II Time 10−2 8.4000e (−003) 5.4393e(−04) 0.21322s 2.6947e(−04) 0.21322s 10−4 1.6621e (−004) 2.6947e(−05) 0.21322s 1.0276e(−06) 0.21322s 10−6 2.7506e (−006) 1.3195e(−08) 0.21321s 6.9756e(−08) 0.21321s Figure 5. Error for Example III 45 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 46 Table 4. Result for Example IV S teps OHBM CASE I CASE II 20 2.943 × 10−5 1.311 × 10−5 1.122 × 10−5 40 1.250×10−6 1.406×10−8 1.465×10−8 80 2.721×10−13 2.545×10−15 1.843×10−15 160 4.262×10−15 2.733×10−18 1.033×10−18 320 7.707×10−18 8.371×10−20 9.572×10−20 640 3.000×10−21 1.721×10−21 1.953×10−21 Computational Time 0.03125s 0.02034s 0.02034s performed better than those of [26]. The comparison of error in Figure 5 clearly shows that, the new methods converge better than the existing method. Example IV: Consider the following linear systems y′1(x) = −21y1 + 19y2 − 20y3 y′2(x) = 19y1 − 21y2 + 20y3 y′3(x) = 40y1 − 40y2 − 40y3  ,  y1 (0) y2 (0) y3 (0)  =  1 0 −1  The exact solution is given as  y1 (x) y2 (x) y3 (x)  =  1 2 ( e−2x + e−40x (cos (40x) + sin (40x)) ) 1 2 ( e−2x − e−40x (cos (40x) + sin (40x)) ) 1 2 e −40x (cos (40x) + sin (40x))  Source: [2] We compared the results of CASE I and CASE II along side OHBM by comparing the maximum relative errors over the three components y1(x), y2(x) and y3(x). As shown in Table 4, the new methods in case I and case II proved to be superior in terms of accuracy and computation time. Figure 6. Error for Example IV Table 5. Result for Example V h [24] ERROR CASE I ERROR CASE II 2π/300 2.58313×10−9 2.03334×10−9 2.03034×10−9 2π/600 1.07326×10−10 2.08320×10−12 2.34020×10−12 2π/1200 3.48507×10−12 5.02513×10−12 5.01013×10−12 46 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 47 Example V: Consider the following problem( y′′(x) = −100y(x) + 99 sin(x) ) , y1 (0) = 1, y ′ (0) = 11, x ∈ [0, 2π] The exact solution is given as y (x) = cos(10x) + sin(10x) + sin(x) Source: [24] Table 6. Result for Example V h [24] Error Case I Error Case 2π/300 2.58313×10−9 2.03334×10−9 2.03034×10−9 2π/600 1.07326×10−10 2.08320×10−12 2.34020×10−12 2π/1200 3.48507×10−12 5.02513×10−12 5.01013×10−12 Results of CASE I and CASE II was compared with [24] as shown in Table 5. The new methods compute favorably with [24] and has better consistency. Figure 7. Error for Example V 4. Conclusion A class of third derivative continuous block methods for the solution of stiff and oscillatory problems is constructed using collocation and interpolation technique. The approximate so- lution adopted for the development of the methods comprises of a combination of polynomial and trigonometric functions. A consistent, zero stable and convergent block methods with con- tinuous coefficients are developed and implemented by writing a code using MATLAB 8.5. The methods are presented in con- tinuous and discrete form. The accuracy of the block methods were tested on some stiff and oscillatory IVP to generate results and compared the results with the results of some existing methods. The results of the new methods compete favorably with the existing methods with better accuracy, consistency and computational time. Thus the methods are consistent, convergent and zero stable. Hence the methods derived are suitable for solving stiff and oscillatory problems and computationally reliable. References [1] A. O. Adesanya, M. R. Odekunle & A. O. Adeyeye, ”Continuous block hybrid predictor corrector method for the solution of y′′ = f (x, y, y′)”, International Journal of Mathematics and Soft Computing 2 (2020) 35. [2] M. O. Ogunniran, Y. Haruna, R. B. Adeniyi & M. O. Olayiwola, ”Op- timized Three-step Hybrid Block Method for Stiff Problems in Ordinary Differential Equations”, Cankaya University Journal of Science and En- gineering 17 (2020) 80. [3] D. Q. Nykamp, ”An introduction to ordi- nary differential equations”, From Math Insight. http://mathinsight.org/ordinary differential equation ntroduction 2020 (2020). [4] D. O. Awoyemi, S. J. Kayode & L. O. Adoghe, ”A six-step continuous multistep method for the solution of general fourth order initial value problems of ordinary differential equations”, Journal of Natural Sciences Research 5 (2015) 2224. [5] S. Adamu, A. A. Danhausa, L. Stephen & B. Williams, ”Two hybrid points block methods for the solution of initial value problems”, J. of NAMP 54 (2020) 7. [6] I., F. Fudziah, Z. A. Sufia, D. J. Yusuf & S. Norazak, ”Block Hybrid Method with Trigonometric-Fitting for Solving Oscillatory Problems”, Sains Malaysiana 47 (2018) 2223. [7] V. O. Marcus & S. Gustaf, ”Classification of Stiffness and Oscillations in Initial Value Problems”, Masters Thesis, Centre for Mathematical Sci- ences, Lund University (2015). [8] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey & D. E. Knuth, ”On the Lambert W Function”, Advances in Computational Mathematics 5 (1996) 329. [9] B. I. Zarina & S. M. Z. Iskandar, ”A Stable Fourth Order Block Backward Differentiation Formulas for Solving Stiff Initial Value Problems”, ASM Sci. J. 6 (2019) 60. [10] B. Bature & H. Musa, ”Stability Analysis of the 2-Point Diagonally Im- plicit Super Class of Block Backward Differentiation Formula with Off- Step Points”, J. Phys Math 8 (2017) 1. [11] N. M. Noor, Z. B. Ibrahim & F. Ismail, ”Numerical solution for stiff ini- tial value problems using 2-point block multistep method”, Journal of Physics: Conf. Series 012017 (2019) 1132. [12] M. O. Udo, R. A. Ademiluyi, A. O. Adesanya & B. O. Ekpenyong, ”An implicit zero stable linear multi-step method for direct solution of gen- eral third order initial value problems of ordinary differential equations”, Journal of Science and Technology Research 6 (2007) 1596. [13] E. A. Areo & R. B. Adeniyi, ”A self-starting linear multistep method for direct solution of initial value problems of second order ordinary differ- ential equations”, International Journal of Pure and Applied Mathematics 8 (2013) 345. [14] A. K. Ezzeddine & G. Hojjati, ”Third derivative multistep methods for stiff systems”, International Journal of Nonlinear Science 14 (2012) 443. [15] J. M. Franco, ”A class of explicit two-step hybrid methods for second- order IVPs”, Journal of Computational and Applied Mathematics 187 (2006) 41. [16] T. A. Anake, D. O. Awoyemi & A. O. Adesanya, ”One-step implicit hy- brid block method for the direct solution of general second order ordinary differential equations”, IAENG International Journal of Applied Mathe- matics 42 (2012). [17] A. O. Adesanya, M. R. Odekunle & M. A. Alkali, ”Three steps block pre- dictor corrector method for the solution of general second order ordinary differential equations”, Engineering Research and Applications 2 (2012) 2297. [18] F. F. Ngwane & S. N. Jator, ”Block hybrid-second derivative method for stiff systems”, International Journal of Pure and Applied Mathematics 80 (2012) 543. [19] R. A. Ademiluyi, M. K. Duromola & B. Bolarinwa, ”Modified block method for the direct solution of initial value problems of fourth order ordinary differential equations”, Australian Journal of Basic and Applied Sciences 8 (2014) 389. [20] K. M. Abualnaja, ”A block procedure with linear multi-step methods us- ing Legendre polynomials for solving ODEs”, Applied Mathematics 6 (2015) 717. [21] F. F. Ngwane & S. N. Jator, ”Solving oscillatory problems using a block hybrid trigonometrically fitted method with two off-step points”, Ninth MSU-UAB Conference on Differential Equations and Computational 47 Kida et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 34–48 48 Simulations. Electronic Journal of Differential Equations Conference 20 (2013) 119. [22] P. L. Ndukum, T. A. Biala, S. N. Jator & R. B. Adeniyi, ”A fourth order trigonometrically fitted method with the block unification implementation approach for oscillatory initial value problems”, International Journal of Pure and Applied Mathematics 103 (2015) 201. [23] F. F. Ngwane & S. N. Jator, ”A family of trigonometrically fitted enright second derivative methods for stiff and oscillatory initial value problems”, Journal of Applied Mathematics 1 (2015). [24] H. Ramos. Z. Kalogiratou, Th. Monovasilis & T. E. Simos, ”A Trigono- metrically Fitted Optimized Two-step Hybrid Block Method for Solving Initial Value Problems of the form y′′ = f (x, y, y′) with Oscillatory Solu- tions”, AIP Conference Proceedings 1 (2015) 1648. [25] H. Musa, M. B. Suleiman & N. Senu, ”Fully Implicit 3-Point Block Extended Backward Differentiation Formula for Stiff Initial Value Prob- lems”, Applied Maths Sciences 6 (2012) 4211. [26] S. A. M. Yatim, Z. B. Ibrahim, K. I. Othman & M. B. Suleiman, ”A Numerical Algorithm for Solving Stiff Ordinary Differential Equations”, Hindawi Publishing Corporation Mathematical Problems in Engineering 11 (2013). [27] H. Musa, M. B. Suleiman, F. Ismail, N. Senu & Z. B. Ibrahim, ”An Im- proved 2-point Block Backward Differentiation Formula for Solving Stiff Initial Value Problems”, AIP Conference Proceedings 211 (2013) 1522. [28] J. Sunday, M. R. Odekunle, A. A. James & A. O. Adesanya, ”Numerical Solution of Stiff and Oscillatory Differential Equations Using a Block In- tegrator”, British Journal of Mathematics and Computer Science 4 (2014) 2471. [29] Z. B. Ibrahim, K. Othman & M. Suleiman, ”Implicit r-point block Back- ward Differentiation Formula for Solving First-order Stiff ODEs”, Ap- plied Mathematics and Computation 186 (2007) 558. [30] A. O. Adesanya, B. Abdulqadri & Y. S. Ibrahim, ”Hybrid one step block method for the solution of third order initial value problems of ordinary differential equations”, International Journal of Applied Mathematics and Computation 6 (2014) 10. [31] R. Abdelrahim & Z. Omar, ”A four-step implicit block method with three generalized off-step points for solving fourth order initial value problem directly”, Journal of King Saud University Science 29 (2017) 401. 48