An Eight Order Two-Step Taylor Series Algorithm for the Numerical Solutions of Initial Value Problems of Second Order Ordinary Differential Equations

Ayodele Olakiitan Owolanke^1*, Ohi Uwaheren², Friday Oghenerukevwe Obarhua³

¹Department of Mathematics, Faculty of Mathematical Sciences, Ondo State University of Science and Technology, Okitipupa, Nigeria.
²Department of Mathematics, Faculty of Physical Sciences, University of Ilorin, Ilorin, Nigeria.
³Mathematical Science Department, Federal University of Technology, Akure, Nigeria.

DOI: 10.4236/oalib.1103486                

Abstract

Our focus is the development and implementation of a new two-step hybrid method for the direct solution of general second order ordinary differential equation. Power series is adopted as the basis function in the development of the method and the arising differential system of equations is collocated at all grid and off-grid points. The resulting equation is interpolated at selected points. We then analyzed the resulting scheme for its basic properties. Numerical examples were taken to illustrate the efficiency of the method. The results obtained converge closely with the exact solutions.

Keywords

Power Series, Collocation and Taylor’s Series Algorithm

1. Introduction

We consider the numerical solution of initial value problem of the form:

$y^{″}=f\left(x,y,y^{\prime }\right);y\left(a\right)=y\left(0\right);y^{\prime }\left(a\right)=\gamma$ (1)

In practice, higher order ordinary differential equations of this form $y^n=f\left(x,y^{\prime },y^{″},\cdots ,y^{n-1}\right)$ , is solved by reducing it to systems of first order differential equation of the form:

$y^{\prime }=f\left(x,y\right),y\left(a\right)=0,f\in c\left[a,b\right],x,y\in {ℝ}^n$ (2)

then an approximate method is applied to solve the resulting Equation (2) as widely discussed by Fatunla [1] and Lambert [2] and Spiegel [3] . The approach does not utilize additional information associated with the specific ordinary differential equation, and consequently, the oscillatory nature of the solution of the differential equation is always neglected. Thus, it would be more efficient to improve on the numerical method so that higher order ordinary differential equations could be solved without having to reduce to systems of first order as suggested by Chakravati and Worland [4] , Dahlquist [5] , sharp and Fine [6] , and Bun and Vasilsyer [7] . Actually, considerable attention has been devoted to solving ordinary differential equation of higher order directly without reduction for instance: methods of linear multistep method (LMM) were considered by Lambert and Watson [8] , Dormand and El-Mikkawy [9] , El-Mikkawy and El- Desouky [10] and Awoyemi [11] [12] [13] [14] . Subsequently, LMM was independently proposed by Kayode [15] , Onumanyi et al. [16] and Adesanya et al. [17] in the predictor-corrector mode, based on collocation method. These authors proposed LMM with continuous coefficients where they adopted Taylor series algorithm to supply the starting values. Also, some notable scholars improve on the predictor-corrector method for solving ordinary differential equations of higher orders, for instance, Jator and Li [18] proposed five-step and four-step methods respectively in which they adopted a continuous LMM to obtain finite difference method. Moreover, Adesanya [19] adopted a method of collocation and interpolation to develop a continuous LMM which is evaluated at different grid points to give discret methods that generate independent solutions. Others that adopted block methods include Badmus and Yahaya [20] . One of the advantages of the method is that it provides direct solution of implicit Linear multistep method without developing separate predictors.

Although some of the aforementioned authors have made use of Taylor series, but little has been said with the use of Taylor series as a major method of implementation. So, Our idea is to use Taylor series algorithm to evaluate

$y_{n+j},{y^{\prime }}_{n+j},j=1,2$ and $y_{n+v},{y^{\prime }}_{n+v},v=\frac{1}{2},\frac{1}{3},\frac{2}{3},\frac{3}{2},\cdots$ and calculate $f^{\prime }$ , $f^{″}$ by the

use partial derivative technique. Thus, two-step hybrid methods in the Taylor series mode are developed to solve second order ordinary differential equations directly.

2. Derivation

In this section, power series is considered as an approximate solution to the general second order problems:

$f\left(x,y,y^{\prime },y^{″}\right)=0;y\left(a\right)=y\left(0\right);y^{\prime }\left(a\right)=\gamma$ (3)

of the form:

$y\left(x\right)={\displaystyle \underset{j=0}{\overset{2k+1}{\sum }}}{a}_j{x}^j$ (4)

The first and second derivative of (3) are respectively given as:

$y^{\prime }\left(x\right)={\displaystyle \underset{j=1}{\overset{2k+1}{\sum }}}j{a}_j{x}^{j-1}$ (5)

$y^{″}={\displaystyle \underset{j=2}{\overset{2k+1}{\sum }}}j\left(j-1\right)a_j{x}^{j-2}$ (6)

Combining (2) and (5), we generate the differential system

${\displaystyle \underset{j=2}{\overset{2k+1}{\sum }}}j\left(j-1\right)a_j{x}^{j-2}=f\left(x,y,y^{\prime }\right),$ (7)

we develop the hybrid scheme using (3) and (5) as interpolation and collocation equations in this work.

Collocating (6) at selected grid and off-grid points, $x=x_{n+1},0\le i\le 2$ and interpolating (3) at selected grid and off-grid points, it results into a system of equations:

${\displaystyle \underset{j=2}{\overset{2k+1}{\sum }}}j\left(j-1\right)a_j{x}^{j-2}=f_{n+i},0\le i\le 2$ (8)

${\displaystyle \underset{j=2}{\overset{2k+1}{\sum }}}{a}_j{x}^j=y_{n+i},0\le i\le 2$ (9)

where, $x_{n+i}=x_n+ih$ , solving Equations ((7) and (8)), ${a^{\prime }}_{j}s$ , yield a method ex-

pressed in the form:

$y_k\left(x\right)={\displaystyle \underset{j=0}{\overset{k}{\sum }}}{\alpha }_j\left(x\right)y_{n+j}+{\displaystyle \underset{j=0}{\overset{k}{\sum }}}{\beta }_j\left(x\right)f_{n+j},$ (10)

where $k=2$ and $f_{n+j}=f\left(x_{n+j},y_{n+j},{y^{\prime }}_{n+j}\right),0\le 2$

It implies

$a_0+a_1{x}_n+a_2{x}_n^2+a_3{x}_n^3+a_4{x}_n^4+a_5{x}_n^5+a_6{x}_n^6+a_7{x}_n^7+a_8{x}_n^8=y_n$ (11)

$a_0+a_1{x}_{n+1}+a_2{x}_{n+1}^2+a_3{x}_{n+1}^3+a_4{x}_{n+1}^4+a_5{x}_{n+1}^5+a_6{x}_{n+1}^6+a_7{x}_{n+1}^7+a_8{x}_{n+1}^8=y_{n+1}$ (12)

$2a_2+6a_3{x}_n+12a_4{x}_n^2+20a_5{x}_n^3+30a_6{x}_n^4+42a_7{x}_n^5+56a_8{x}_n^6=f_n$ (13)

$2a_2+6a_3{x}_{n+\frac{1}{3}}+12a_4{x}_{n+\frac{1}{3}}^2+20a_5{x}_{n+\frac{1}{3}}^3+30a_6{x}_{n+\frac{1}{3}}^4+42a_7{x}_{n+\frac{1}{3}}^5+56a_8{x}_{n+\frac{1}{3}}^6=f_{n+\frac{1}{3}}$ (14)

$2a_2+6a_3{x}_{n+\frac{2}{3}}+12a_4{x}_{n+\frac{2}{3}}^2+20a_5{x}_{n+\frac{2}{3}}^3+30a_6{x}_{n+\frac{2}{3}}^4+42a_7{x}_{n+\frac{2}{3}}^5+56a_8{x}_{n+\frac{2}{3}}^6=f_{n+\frac{2}{3}}$ (15)

$2a_2+6a_3{x}_{n+1}+12a_4{x}_{n+1}^2+20a_5{x}_{n+1}^3+30a_6{x}_{n+1}^4+42a_7{x}_{n+1}^5+56a_8{x}_{n+1}^6=f_{n+1}$ (16)

$2a_2+6a_3{x}_{n+\frac{4}{3}}+12a_4{x}_{n+\frac{4}{3}}^2+20a_5{x}_{n+\frac{4}{3}}^3+30a_6{x}_{n+\frac{4}{3}}^4+42a_7{x}_{n+\frac{4}{3}}^5+56a_8{x}_{n+\frac{4}{3}}^6=f_{n+\frac{4}{3}}$ (17)

$2a_2+6a_3{x}_{n+\frac{5}{3}}+12a_4{x}_{n+\frac{5}{3}}^2+20a_5{x}_{n+\frac{5}{3}}^3+30a_6{x}_{n+\frac{5}{3}}^4+42a_7{x}_{n+\frac{5}{3}}^5+56a_8{x}_{n+\frac{5}{3}}^6=f_{n+\frac{5}{3}}$ (18)

$2a_2+6a_3{x}_{n+2}+12a_4{x}_{n+2}^2+20a_5{x}_{n+2}^3+30a_6{x}_{n+2}^4+42a_7{x}_{n+2}^5+56a_8{x}_{n+2}^5=f_{n+2}$ (19)

Writing these system of equations in matrix form:

$\left(\begin{array}{ccccccccc}1& x_n& x_n^2& x_n^3& x_n^4& x_n^5& x_n^6& x_n^7& x_n^8\\ 1& x_{n+1}& x_{n+1}^2& x_{n+1}^3& x_{n+1}^4& x_{n+1}^5& x_{n+1}^6& x_{n+1}^7& x_{n+1}^8\\ 0& 0& 2& 6x_n& 12x_n^2& 20x_n^3& 30x_n^4& 42x_n^5& 56x_n^6\\ 0& 0& 2& 6x_{n+1/3}& 12x_{n+1/3}^2& 20x_{n+1/3}^3& 30x_{n+1/3}^4& 42x_{n+1/3}^5& 56x_{n+1/3}^6\\ 0& 0& 2& 6x_{n+2/3}& 12x_{n+2/3}^2& 20x_{n+2/3}^3& 30x_{n+2/3}^4& 42x_{n+2/3}^5& 56x_{n+2/3}^6\\ 0& 0& 2& 6x_{n+1}& 12x_{n+1}^2& 20x_{n+1}^3& 30x_{n+1}^4& 42x_{n+1}^5& 56x_{n+1}^6\\ 0& 0& 2& 6x_{n+4/3}& 12x_{n+4/3}^2& 20x_{n+4/3}^3& 30x_{n+4/3}^4& 42x_{n+4/3}^5& 56x_{n+4/3}^6\\ 0& 0& 2& 6x_{n+5/3}& 12x_{n+5/3}^2& 20x_{n+5/3}^3& 30x_{n+5/3}^4& 42x_{n+5/3}^5& 56x_{n+5/3}^6\\ 0& 0& 2& 6x_{n+2}& 12x_{n+2}^2& 20x_{n+2}^3& 30x_{n+2}^4& 42x_{n+2}^5& 56x_{n+2}^6\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\\ a_7\\ a_8\end{array}\right)=\left(\begin{array}{c}{y}_n\\ y_{n+1}\\ f_n\\ f_{n+1/3}\\ f_{n+2/3}\\ f_{n+1}\\ f_{n+4/3}\\ f_{n+5/3}\\ f_{n+2}\end{array}\right)$ (20)

Using Gaussian elimination method, the unknown coefficients ${a^{\prime }}_{j}s$ can be obtained. Putting ${a^{\prime }}_{j}s$ back into (3) gives (10):

The coefficients ${\alpha }_{i^{\prime }s}\left(t\right)$ , ${\beta }_{j^{\prime }s}\left(t\right)$ are continuous coefficients obtained using the transformation $t=\frac{1}{h}\left(x-x_{n+k-1}\right)$ , $t\in \left(\mathrm{0,1}\right]$

$\frac{\text{d}t}{\text{d}x}=\frac{1}{h}.$

Then simplifying the continuous ${\alpha }_{j^{\prime }s}$ , ${\beta }_{j^{\prime }s}$ , and taking their first derivatives, we have:

${\alpha }_0{\left(t\right)}^{\prime }=-\frac{1}{h}$ ,

${\alpha }_1{\left(t\right)}^{\prime }=-\frac{1}{h}$ ,

${\beta }_0{\left(t\right)}^{\prime }=\frac{47h}{13440}$ ,

${\beta }_{\frac{1}{3}}{\left(t\right)}^{\prime }=\frac{327h}{2240}$ ,

${\beta }_{\frac{2}{3}}{\left(t\right)}^{\prime }=\frac{111h}{890}$ ,

${\beta }_1{\left(t\right)}^{\prime }=\frac{1088h}{3360}$ ,

${\beta }_{\frac{4}{3}}{\left(t\right)}^{\prime }=\frac{93h}{640}$ ,

${\beta }_{\frac{5}{3}}{\left(t\right)}^{\prime }=\frac{1095h}{2240}$ ,

${\beta }_2{\left(t\right)}^{\prime }=\frac{1359h}{13440}$ .

Then, putting $t=1$ gives:

$\begin{array}{c}{y}_{n+2}=2y_{n+1}+y_n+\frac{h^2}{6720}\{47f_{n+2}+810f_{n+\frac{5}{3}}+1377f_{n+\frac{4}{3}}\\ +2252f_{n+1}+1377f_{n+\frac{2}{3}}+810f_{n+\frac{1}{3}}+857f_n\}\end{array}$ (21)

its first derivative

$\begin{array}{c}{y^{\prime }}_{n+2}=\frac{1}{h}\left[y_{n+1}-y_n\right]+\frac{h^2}{6720}\{1359f_{n+2}+6570f_{n+\frac{5}{3}}+1953f_{n+\frac{4}{3}}\\ +4352f_{n+1}+1665f_{n+\frac{2}{3}}+1962f_{n+\frac{1}{3}}+47f_n\}\end{array}$ (22)

with the order $p=8$ , error constant $C_{10}=-0.0069941$ , and interval of absolute stability $X\left(\Theta \right)=\left(-14.1608,0\right)$ Implementation of the method using Taylor series algorithm to evaluate

$y_{n+j},{y^{\prime }}_{n+j},y_{n+v},{y^{\prime }}_{n+v},f_{n+v},f_{n+j}$ ,

where, $j^{\prime }s=1,2$ and $v^{\prime }s=\frac{1}{3},\frac{2}{3},\frac{4}{3},\frac{5}{3}$ and,

$f_{n+v}=f\left(x_{n+v},y_{n+v},{y^{\prime }}_{n+v}\right),$

such that

$y_{n+v}=y_n+vh{y^{\prime }}_n+\frac{{\left(vh\right)}^2}{2!}{f}_n+\frac{{\left(vh\right)}^3}{3!}{f^{\prime }}_n+\frac{{\left(vh\right)}^4}{4!}{f^{″}}_n+\cdots$ (23)

and,

${y^{\prime }}_{n+v}={y^{\prime }}_n+vh{f}_n+\frac{{\left(vh\right)}^2}{2!}{f^{\prime }}_n+\frac{{\left(vh\right)}^3}{3!}{f^{″}}_n+\frac{{\left(vh\right)}^4}{4!}{f^{‴}}_n+\cdots$ (24)

Also,

$f_{n+j}=y^{″}\left(x_n+jh\right)=f_n+jh{f^{\prime }}_n+\frac{{\left(jh\right)}^2}{2!}{f^{″}}_n+\cdots$ (25)

$f_n=f\left(x_n,y_n,{y^{\prime }}_n\right)$

$f^{\left(i\right)}=f^{\left(i\right)}\left(x_n,y_n,{y^{\prime }}_n\right),i=1,2$

Finding the partial derivative $f^{\prime }\mathrm{,}{f}^{″}\mathrm{,}\cdots$ as follows

$\frac{\text{d}f}{\text{d}x}=f^{\prime }=\frac{\delta f}{\delta x}+\frac{\delta f}{\delta x}{y}^{\prime }+\frac{\delta f}{\delta y^{\prime }}f$ (26)

$f^{″}=\frac{{\text{d}}^{2}f}{\text{d}{x}^2}=2\left(A{y}^{\prime }+Bf\right)+Cf{y}^{\prime }+D+E,$ (27)

where,

$A=\frac{{\partial }^{2}f}{\partial x\partial y}+f\frac{{\partial }^{2}f}{\partial y\partial y^{\prime }}$ (28)

$B=\frac{{\partial }^{2}f}{\partial x\partial y^{\prime }}$ (29)

$C=\frac{\partial f}{\partial x}+y^{\prime }\frac{\partial f}{\partial y}+f\frac{\partial f}{\partial y^{\prime }}$ (30)

$D=\frac{{\partial }^{2}f}{\partial x^2}+{\left(y^{\prime }\right)}^2\frac{{\partial }^{2}f}{\partial y^2}+f^2\frac{{\partial }^{2}f}{\partial {\left(y^{\prime }\right)}^2}$ (31)

$E=f\frac{\partial f}{\partial y}$ (32)

2.1. Analysis of the Properties of the Scheme

We shall consider the analysis of the basic properties of our methods which includes the order, the region of absolute stability and the zero stability of the methods.

2.2. Order of Accuracy of the Method

The local truncation error with k-step linear multistep m method which is in line with Lambert (1973), is taken to be linear difference operator $\mathcal{l}$ defined by

$\mathcal{l}\left[y\left(x\right);h\right]={\displaystyle \underset{j=0}{\overset{k}{\sum }}}\left[{\alpha }_{j}y\left(x_n+j\right)-h{\beta }_{j}y\left(x_n+j\right)\right]$ (33)

Thus, expanding (21) as Taylor series about point x and comparing coefficients of $h^k$ , the scheme will be of order $p=8$ with error constant $C_{p+2}=-0.0069941$

$L\left[y\left(x\right),h\right]=C_{0}y\left(x_n\right)+C_1{y}^{\prime }\left(x_n\right)+C_2{y}^{″}\left(x_n\right)+\cdots +C_p{y}^p\left(x_n\right),$ (34)

where $C_p,p=0,1,\cdots$ are the constant coefficients given as:

$\begin{array}{l}{C}_0={\displaystyle \underset{j=0}{\overset{k}{\sum }}}{\alpha }_j\\ C_1={\displaystyle \underset{j=0}{\overset{k}{\sum }}}j{\alpha }_j\text{and}\\ C_p=\frac{1}{p!}\left[{\displaystyle \underset{j=0}{\overset{k}{\sum }}}j{\alpha }_j-p\left(p-1\right)\left({\displaystyle \underset{j=0}{\overset{k}{\sum }}}{j}^{p-1}{\beta }_j+{\displaystyle \underset{j=0}{\overset{k}{\sum }}}{q}^{p-1}{\beta }_{qj}\right)\right]\end{array}\}$ (35)

In line with [2] , k-step, linear multistep (21) has order p if $C_0=C_1=\cdots =C_{p-1}=C_p$ and $C_{p+1}\ne 0$ , where, $C_{p+1}\ne 0$ is the error constant. Subjecting our schemes to equations 35, it is therefore established that linear multistep scheme is of order $p=8$ , relatively small error constant −0.0069941.

2.3. Consistency of the Scheme

A linear multistep method is consistent if the following conditions are satisfied:

1) The order $p\ge 1$ .

2) $p\left(1\right)=0,p^{\prime }\left(1\right)=\sigma \left(1\right)$ .

3) ${\displaystyle {\sum }_{j=0}^k}{\alpha }_j=0$ .

4) ${\displaystyle {\sum }_{j=0}^k}j{\alpha }_j={\displaystyle {\sum }_{j=0}^k}{\beta }_j$ .

2.4. Zero Stability of the Method

Equation (21) has its first characteristic polynomial to be:

$\rho \left(r\right)=r^2-2r+1$ (36)

The method is zero stable since they have roots $r=1$ twice.

2.5. Region of Absolute Stability of the Method

In order to establish the region of absolute stability, we apply the boundary locus method as in [2] . The method implies that

$\theta =\frac{\rho \left(r\right)}{\delta (\; r\; )}$

where,

$r={\text{e}}^{i\theta }=\cos \left(\theta \right)+i\sin (\; \theta \; )$

From scheme (21), we have: $\rho \left(r\right)=r^2-2r+1$

and

$\sigma \left(r\right)=\frac{1}{6720}\left[47r^2+810r^{\frac{5}{3}}+1377r^{\frac{4}{3}}+2252r+1377r^{\frac{2}{3}}+810r^{\frac{1}{3}}+857\right]$

so that

$h\left(\theta \right)=\frac{\rho \left({\text{e}}^{i\theta }\right)}{\delta \left({\text{e}}^{i\theta }\right)}$

which implies

$h\left(\theta \right)=\frac{1}{6720}\left[47r^2+810r^{\frac{5}{3}}+1377r^{\frac{4}{3}}+2252r+1377r^{\frac{2}{3}}+810r^{\frac{1}{3}}+857\right]$ (37)

$\begin{array}{c}h\left(\theta \right)=\left[\cos \left(2\theta \right)+i\sin \left(2\theta \right)-2\cos \left(\theta \right)-2i\sin \left(\theta \right)+1\right]\\ \times 6720(47\cos \left(2\theta \right)+47i\sin \left(2\theta \right)+810\cos \left(\frac{5\theta }{3}\right)+810i\sin \left(\frac{5\theta }{3}\right)\\ +1377\cos \left(\frac{4\theta }{3}\right)+1377i\sin \left(\frac{4\theta }{3}\right)+2252\cos \left(\theta \right)+2252i\sin \left(\theta \right)\\ +1377\cos \left(\frac{2\theta }{3}\right)+{1377i\sin \left(\frac{2\theta }{3}\right)+810\cos \left(\frac{2\theta }{3}\right)+810i\sin \left(\frac{2\theta }{3}\right)+47)}^{-1}\end{array}$

Considering the values of $\theta$ for $0\le \theta \le 180$ at intervals of ${30}^{\theta }$ gives the region of absolute stability to be $\left(-\mathrm{14.1608,0}\right).$

3. Numerical Experiments

We test the accuracy of the proposed scheme on some numerical problems, and the results are compared with existing methods.

Problem 1:

$y^{″}=x{\left(y^{\prime }\right)}^2,y\left(0\right)=1,y^{\prime }\left(0\right)=0.5,h=\frac{0.1}{32}$ (38)

Exact solution

$y\left(x\right)=1+\frac{1}{2}{\log }_{10}\left(\frac{2+x}{2-x}\right)$

The numerical results of the problem is shown in Table 1, and is compared with Awoyemi and kayode (2005) of order 8 in Table 2.

Problem 2:

$y^{″}=\left(-6/x\right)y^{\prime }-\left(4/{\left(x\right)}^2\right)y,y\left(1\right)=1,y^{\prime }\left(1\right)=1,h=120$ (39)

Table 1. Results and errors for problem (1).

Table 2. Results and errors for problem (2).

Note: YEX = Yexact, YC = Ycomputed, ERRNew = Error in new method.

Exact solution

$y\left(x\right)=1-{\text{e}}^x$

The numerical results of the problem is shown in Table 2.

4. Conclusion

A Linear Multistep method which implements a Taylor’s series algorithm is developed for the direct solution of general second order initial value problems of

ordinary differential equations without reduction to systems of first order differential equation. The derivatives of continuous scheme to any order were computed implementing Taylor’s series algorithm. The accuracy of the method was tested with two test problems, and results were compared with Awoyemi and Kayode [11] of order (8).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Fatunla, S.O. (1988) Numerical Methods for Initial Value Problems for Ordinary Differential Equations. Academy Press, New York, 295.
[2]	Lambert J.D. (1973) Computational Methods in Ordinary Differential Equations. John Wiley and Sons, New York, 278.
[3]	Murray, S.R. (1971) Theory and Problems of Advance Mathematics for Engineers and Scientist. Mc Graw Hill, Inc., New York.
[4]	Chakravati, P.C. and Worland, P.B. (1971) A Class of Self Starting Methods for the Numerical Solution of Ordinary Differential Equation. BIT, 11, 368-383. https://doi.org/10.1007/BF01939405
[5]	Dahlquist, G. (1959) Convergence and Stability in the Numerical Integration of Ordinary Differential Equation. Mathematics Scandinavia, 4, 33-53. https://doi.org/10.7146/math.scand.a-10454
[6]	Sharp and Fine (1992) Some Nystrom Pairs for the General Second Order Value Problem. Journal of Computation and Applied Mathematics, 42, 279-291. https://doi.org/10.1016/0377-0427(92)90081-8
[7]	Bun, R.A. and Vasilsyer, Y.D. (1990) A Numerical Methods for Solving Differential Equations of Any Order. Computational Mathematics and Mathematical Physics, 32, 317-330.
[8]	Lambert, J.D. and Wastson, A. (1976) Symmetric Multistep Method for Periodic Initial Value Problem. Journal of the Institute of Mathematics Application, 18, 189- 202. https://doi.org/10.1093/imamat/18.2.189
[9]	Dormand, J.R., El-Mikkawy, M.E.A. and Prince, P.J. (2003) Families of Runge Kutta Nystrom Formula, IMA. Journal of Numerical Analysis, 7, 235-250. https://doi.org/10.1093/imanum/7.2.235
[10]	El-Mikkawy, M.E. and El-Desouky, R. (2003) A New Optimized Non FSAL Embedded Runge-Kutta-Nystrom Algorithms of Order 6 and 4 in Six Stages. Applied Mathematics and Computation, 145, 33-43. https://doi.org/10.1016/S0096-3003(02)00436-8
[11]	Awoyemi, D.O. and Kayode, S.J. (2005) A Maxima Order Collocation Method for Initial Value Problems of General Second Order Ordinary Differential Equation. Proceeding of the Conference Organized by the National Mathematical Center, Abuja.
[12]	Awoyemi, D.O. (1999) A Class of Continuous Method for General Second Order Initial Value Problems in Ordinary Differential Equations. International Journal of Computer Mathematics, 72, 29-37. https://doi.org/10.1080/00207169908804832
[13]	Awoyemi, D.O. (2001) A New Sixth Order Algorithms for General Second Order Ordinary Differential Equation. International Journal of Computer Mathematics, 77, 117-124. https://doi.org/10.1080/00207160108805054
[14]	Awoyemi, D.O. (2005) An Algorithm Collocation Methods for Direct Solution of Special and General Fourth Order Initial Value Problems of Ordinary Differential Equation. Journal of Nigerian Association of Mathematical Physics, 6, 271-238.
[15]	Kayode, S.J. (2009) An Order Zero Stable Method for Direct Solution of Fourth Order Ordinary Differential Equation. American Journal of Applied Sciences, 5, 1461-1466.
[16]	Onumanyi, P., Sirisons, U.W. and Dauda, Y. (2001) Toward Uniformly Accurate Continuous Finite Difference Approximation of Ordinary Differential Equation. Bayero Journal of Pure and Applied Sciences, 1, 5-8.
[17]	Adesanya, A.O., Anake, T.A. and Udoh, M.O. (2008) Improved Continuous Method for Direct Solution of General Second Order Ordinary Differential Equation. Journal of Nigeria Association of Mathematical Physics, 13, 59-62.
[18]	Jator, S.N. and Li, J. (2009) A Self Stationary Linear Multistep Method for a Direct Solution of General Second Order Initial Value Problem. International Journal of Computer Mathematics, 86, 817-1165.
[19]	Adesanya, A.O. (2011) Block Methods for Direct Solution of General Higher Order Initial Value Problem of Ordinary Differential Equations. PhD Thesis, Federal University of Technology, Akure.
[20]	Badmus, A.M. and Yahaya, Y.A. (2009) An Accurate Uniform Order 6 Block Method for Direct Solution of General Second Order Ordinary Differential Equations. The Pacific Journal of Science and Technology, 10, 273-277.