Block Methods for Solving Second Order Ordinary Differential Equation Using Shifted Linear Multistep Formulas
- Utalor I.K
- Durojaye M.O
- Ajie Jame I.
- 218-229
- Aug 29, 2025
- Mathematics
Block Methods for Solving Second Order Ordinary Differential Equation Using Shifted Linear Multistep Formulas
Utalor I.K*, Durojaye M.O, Ajie Jame I.
Department of Mathematics science, University of Abuja. FCT
*Corresponding Author
DOI: https://doi.org/10.51584/IJRIAS.2025.100800019
Received: 25 July 2025; Accepted: 29 July 2025; Published: 29 August 2025
ABSTRACT
In this paper, a new block numerical integrator with characteristics of linear multi-step formulas (LMFs) is derived, analyzed, and numerically applied to solve second-order ordinary differential equations (ODEs). It was done by applying a shift operator to two linear multi-step formulae with low Local Truncation Error (LTE). The mathematical derivation of the Linear Multi-step Formulas (LMF) used in the method is based on the method of undetermined coefficients. The analysis of the basic properties of the method shows that the scheme is consistent, convergent, and zero stable. The Numerical experimentation and comparative analysis with existing methods show that our scheme is efficient.
Keywords: One-block Numerical Integrator Method; Shift Operator, Undetermined Coefficients, General Second Order, Initial Value Problems, Self-starting.
INTRODUCTION
Linear multistep methods constitute a powerful class of numerical procedures for solving a second order equation of the form
y” (t)=f(t, y(t), y’ (t)), y (t0) = y0, y’ (t0)=y1, t∈[a,b] (1)
Many Engineers and Scientists spend much of their time working with the models of the world around them. Usually, these mathematical models are developed to help in the understanding of physical phenomena. These models often yield equations that contain some derivatives of an unknown function of one or several variables. Such equations are called differential equations. Differential equations do not only arise in the physical sciences but also in diverse fields such as economics, medicine, and even in areas of biology, just to mention a few. Most often, many of these problems are difficult to be solved by analytical approach. As a result of this, approximate numerical integration methods are routinely used. Conventionally, higher order ODEs are solved by reducing it to a system of first order equations (Lambert, 1973; Fatunla,1988; Brugnano and Trigiante, 1998). The approach results to more cumbersome computation which leads to more errors during the integration process (Adesanya et al. 2009; Vigo-Aguilar et al. 2006).Some attempts have been made to solve problem (1) directly without reduction to a first order systems of equations (Adesanya etal., 2008; Badmus & Yahaya, 2009). Authors such as Okunnuga (2008),Onumanyi et al. (1999), and Omar & Kuboye (2015) to mention a few, have presented continuous linear multistep techniques. To obtain beginning values for their approaches, these authors used the predictor-corrector, block methodology and Taylor series expansion. According to Adesanya, (2011), the predictor corrector method is expensive because subroutines are difficult to construct due to the specific methods required to give starting values and change the step size, resulting in longer computer time and more human labor. The predictors and corrector are not in the same order and as a result of this, the method accuracy is poor. A number of authors including Badmus,et al.(2009), Omar,et al.(2015), Onumanyi et al. (1999) and some others have developed the block technique. The undetermined coefficients method was used to construct the coefficients of Boundary Value Methods (BVM) , Linear Multi-step Formula (LMF) with aid of shift operator that possess good stability properties suitable for efficient solutions of stiff initial value problems. These methods have several advantages among which are, the fact that, they do not require any starting value to kick-start their implementation and they are self-starting. It was implemented on first order ordinary differential equations, (Ajie et al.2019; Ajie et al., 2020). This paper is an extension of the work by constructing a method from the perspective of undetermined coefficients method using two Linear Multi-step Formula (LMF) of low Local Truncation Error (LTE ) out of the main method and its derivative with the aid of shift operator. We aimed to obtain a block method with the lowest LTE to formulate it and to apply it to solve second order ordinary differential equations directly without reducing it to first order ordinary differential equations.
The present work is outlined as follows: In Section 2, we present the method for solving second order ODEs. The characteristics of the developed formulas are analyzed in Section 3. In Section 4 and 5, we present the numerical results of some test problems to show the efficiency and reliability of the proposed technique. Conclusion is outlined in Section 6.
DERIVATION OF THE SCHEMES
Consider an approximate solution to the general second order ODEs initial value problems by a self-starting block method. The continuous coefficients
\(\{\alpha_j(t)\}_{j=0}^k, \{\beta_i(t)\}_{i=0}^k, \{\gamma_i(t)\}_{i=0}^k\) for the composing LMF are determined by imposing order condition on linear multi-step formula (LMF) and using the method of undetermined coefficients.
Proposition
On the self-starting block methods, let \( l, m (\geq 2) \) be integers, also let \( m \) denote the number of \( k \)-step LMF in a composite scheme having order of at least \( P_{(k-j)} (\geq 1), j,k(\geq 1)\). Then, the technique of deriving block methods for solving a second order problem is given by applying shift-operator \( l \) times, where \( l \) is given as:
\[
l = \frac{2k-m}{m-2} \quad \text{(Utalor et al., 2025).}
\]
PROOF:
The E-operator is effectively applied \( l \) times on the system of LMF. Thus, there are \( 2(k+l) \) unknown solution points captured in the block of solution. By this block definition,
\[
A_1 h^\lambda Y_m^{(n)} = h^\lambda \sum_{i=0}^k \left[ A_0 Y_{m-i} + h^\mu \left(\sum_{i=1}^k B_1 F_m + B_0 F_{m-i}\right)\right],
\quad \det(A_1) \neq 0
\]
is realized if the coefficients \( A_1, A_0, B_1, B_0 \) are square matrices of dimensions \( 2(k+l) \times 2(k+l) \) for a fixed \( m \), and \( Y_m, Y_{m-i}, F_m, F_{m-i}; m=0,1,2,\dots \) are vectors as specified below:
\[
h^\lambda Y_m^{(n)} = (y_{n+1}, y_{n+2}, y_{n+3}, \dots, h y’_{n+1}, h y’_{n+2}, h y’_{n+3}, \dots, h^2 y^{”}_{n+1}, h^2 y^{”}_{n+2}, \dots, h^2 y^{”}_{n+m})^T
\]
\[
h^\lambda Y_{m-i} = (y_{n+k-1}, y_{n+k-2}, y_{n+k-3}, \dots, y_n,
h y’_{n+k-1}, h y’_{n+k-2}, h y’_{n+k-3}, \dots, h y’_n,
h^2 y^{”}_{n+k-1}, h^2 y^{”}_{n+k-2}, \dots, h^2 y^{”}_n)^T
\]
\[
F_{m-i} = (f_{n+k-1}, f_{n+k-2}, f_{n+k-3}, \dots, f_n, \dots, G_{n+k-1}, G_{n+k-2}, G_{n+k-3}, \dots, G_n, \dots)^T
\]
\[
F_m = (f_{n+1}, f_{n+2}, f_{n+3}, \dots, f_{n+k}, G_{n+1}, G_{n+2}, G_{n+3}, \dots, G_{n+k}) \tag{2}
\]
This simply implies that:
\[
m + ml = 2(k+l)
\]
when \(k\) is chosen such that \(l\) is an integer given as:
\[
l = \frac{2k-m}{m-2}, \quad m \geq 4, \; k \geq 2, \; l \geq 0 \tag{3}
\]
Where
(4)
In this case, the number of equations is m=4 . (Utalor et al., 2025).
Main Formulas (k=3)
Let us consider the Linear Multi-step formula (LMF) of the form
Here, the number of unknowns \( y_{(n+j)}, y’_{(n+j)}, \; j=1,2,3 \) is greater than the number of equations in (9). Originally, there are four equations in equation (9).
By applying the theory given in equation (3) on equation (9), where \( k=3, m=4 \), the result obtained requires the shift operator once. Shifting once, additional four equations are added which makes it eight equations now with eight unknowns:
\[
\{ y_{(n+j)}, y’_{(n+j)} \}, \quad j=1,2,3,4
\]
Since the number of unknowns and the number of equations are equal, then the coefficients of the resultant block method (4) after the shift operator application in vector form are given below:
FUNDAMENTAL PROPERTIES OF THE METHOD
Order and Local Truncation Error
The linear differential operator \( L[z(x); h] \) associated with the block (10) is defined by:
\[
L[y(t);h]=A_1 h^\lambda Y_m^{(n)}-h^\lambda \sum_{i=0}^k \left(A_0 Y_{m-i}+h^\mu \left(\sum_{i=1}^k B_1 F_m+B_0 F_{m-i}\right)\right) \tag{11}
\]
Expanding (10) using Taylor series, we obtain
\[
L[y(t);h]=C_0 y(t)+C_1 h y'(t)+C_2 h^2 y”(t)+\dots+C_q h^q y^{(q)}(t)+\dots \tag{12}
\]
where \( C_q, q=0,1,2,\dots \) are constants given in terms of \(\alpha_j, \beta_j\) and \(\lambda_j\). So that
\[
L[y(t);h]=C_{p+2} h^{p+2} y^{(p+2)}(t)+O(h^{p+3}), \quad p=8
\]
with
\[
C_0=C_1=C_2=\dots=C_{p+1}=0, \quad C_{10}\neq 0=C_{p+2}.
\]
In this case, \(p\) is the order and \(C_{p+2}\) is the error constant (Lambert, 1973).
The error constants \( C_{p+2}^{(T_j^{(\nu)})} \) for \( j=1,2 \) and \( \nu=0,1 \) are given below:
\[
C_{p+2}^{(T_j^{(\nu)})}=\left(-\frac{1}{14025},\;-\frac{137}{6283200},\;-\frac{2}{32725},\;-\frac{163}{4188800}\right)^T
\]
Consistency
Definition 1: A block method is said to be consistent if its order is greater than one (Lambert, 1973). It follows that for all the formulae, the order is \(p=8\). Since the order of the formulae is greater than one, it means that they are consistent.
Zero Stability
Definition 2: The implicit block method (10) is said to be zero stable if the roots \(z_s, s=1,\dots,n\) of the first characteristic polynomial \(\bar{\rho}(z)\), defined by
\[
\bar{\rho}(z)=\det[z\bar{A}_1 – A_0]
\]
satisfy \(|z_s|\leq 1\) and every root with \(|z_s|=1\) has multiplicity not exceeding two in the limit as \(h \to 0\) (Henrici, 1962).
Using the definitions, the method in (10) may be rewritten in a more appropriate vector form to study zero-stability as:
\[
A_1 y_m^{(n)} – A_0 y_{m-1} = 0
\]
where
\[
y_m^{(n)} = (y_{n+1}, y_{n+2},\dots,y_{n+k})^T, \quad y_{m-i} = (y_{n-1}, y_{n-2},\dots,y_n)^T
\]
and \(A_1, A_0\) are constant matrices given by
\[
\ddot{\rho}(z) = \det\!\Bigg[z
\begin{bmatrix}
0 & 1 & 0 & 0 \\
1 & 0 & 0 & 0 \\
1 & 0 & 1 & 0 \\
1 & 1 & 0 & 0
\end{bmatrix}
–
\begin{bmatrix}
0 & 0 & 0 & 1 \\
0 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix}\Bigg]
\]
\[
=-2z^3=0 \quad \therefore z=0
\]
Convergence
Theorem: Consistency and zero stability are sufficient conditions for a linear multistep method to be convergent (Henrici, 1962). Since the method (10) is consistent and zero stable, it implies that the method is convergent.
NUMERICAL ILLUSTRATIONS
In order to ascertain the efficiency of the developed method, numerical experiments of some problems are performed and the results compared with earlier ones in the literature. The accuracy is measured by using the formula:
\[
E_{rc} = \|y(t_j) – y_j\|
\]
where \(E_{rc}\) denotes the absolute error at the considered node, \(y(t_j)\) and \(y_j\) are exact and approximate solutions of the problems, respectively. The computed results for the three problems using the methods proposed are presented in tables and graphically.
Problem 1: Application in Physical Sciences and Engineering – Mass Spring Motion
Consider an application where a 128 lb weight is attached to a spring having a spring constant of 64 lb/ft. The weight is set in motion with no initial velocity by displacing it 6 inches above the equilibrium position and by simultaneously applying to the weight an external force:
\[
F(t) = 8\sin(4t)
\]
Assuming no air resistance, compute the subsequent motion of the weight at \(0.01 \leq t \leq 0.10\).
Modeling this problem into a mathematical equation and applying the proposed method to compute the motion of the weight attached to the spring. The following parameters were considered:
\(m=128,\;k=64,\;b=0,\;F(t)=0.\)
Thus, the application problem is written mathematically as:
\[
\frac{d^2x}{dt^2}+16x = 2\sin(4t), \quad x(0)=-\tfrac{1}{2}, \quad x'(0)=0
\]
with the analytical solution given below:
\[
x(t) = -\tfrac{1}{2}\cos(4t)+\tfrac{1}{16}\sin(4t)-\tfrac{t}{4}\cos(4t) \quad \text{(Ogunware et al., 2023).}
\]
The new method is applied to solve problem 1 and then compared the result with Ogunware et al., (2023) and Skwame et al., (2018) as shown below:
Table 1. Comparison of Absolute errors for problem 1( h = 0.1)
t | Exact solution | Computed solution | Error | Ogunware et al., (2023), | Skwame et al., (2018), |
0.1 | -0.499598720210477 | -0.499598720223385 | 1.2908e-11 | 3.0489e-11 | 1.6621e-09 |
0.2 | -0.498390193309750 | -0.498390193339769 | 3.0019e-11 | 6.1223e−11 | 1.1587e-08 |
0.3 | -0.496368369740280 | -0.496368369787363 | 4.7083e-11 | 9.2137e-11 | 2.9743e-08 |
0.4 | -0.493528526608179 | -0.493528526672251 | 6.4072e-11 | 1.2316e-10 | 8.6076e-08 |
0.5 | -0.489867287968945 | -0.489867288029452 | 6.0507e-11 | 1.5423e-10 | 9.0505e-08 |
0.6 | -0.485382642897099 | -0.485382642947256 | 5.0157e-11 | 1.8526e-10 | 1.3291e-07 |
0.7 | -0.480073961290567 | -0.480073961330294 | 3.9727e-11 | 2.1620e-10 | 1.8317e-07 |
0.8 | -0.473942007364362 | -0.473942007393595 | 2.9233e-11 | 2.4698e-10 | 2.4110e-07 |
0.9 | -0.466988950792028 | -0.466988950739920 | 5.2108e-11 | 2,7751e-10 | 3.0654e-07 |
1 | -0.459218375457224 | 0.459218375300726 | 1.56498e-10 | 3.0773e-10 | 3.7922e-07 |
Fig. 1. Plots of analytical and computed solution for Problem (1).
It shows good agreement between the numerical and exact solutions
Problem 2. Consider a highly stiff initial value problem
\[
f(t, y, y’) = 100y,\qquad y(0)=0,\qquad y'(0)=-10
\]
Exact solution
\[
y(t) = e^{-10t}, \qquad \text{with } h=\frac{1}{100}
\]
as described in(Opeyemi et al, 2024; Akeremale et al, 2024).
The new method is applied to solve problem 2 and then compared the result with
(Opeyemi et al, 2024; Akeremale et al, 2024) as shown below:
Table 2. Comparing the absolute errors in the new methods to others(problem 2)
t | Exact solution | Computed solution | Error | Opeyemi et al., (2024), | Akeremale et al., (2024), |
0.01 | 0.904837418035960 | 0.904837418035940 | 2.00000e-14 | 1.791916e-9 | 2.04 e-7 |
0.02 | 0.818730753077982 | 0.818730753077918 | 6.40000e-14 | 4.313791e-9 | 1.85e-7 |
0.03 | 0.740818220681718 | 0.740818220681572 | 1.46000e-13 | 6.867877e-9 | 1.68e-7 |
0.04 | 0.670320046035639 | 0.670320040300624 | 5.735015e-9 | 8.831792e-9 | 1.52 e-7 |
0.05 | 0.606530659712633 | 0.60653064316728 | 1.654535e-8 | 2.292857e-8 | 1.38e-7 |
0.06 | 0.548811636094027 | 0.54881162053084 | 1.556317e-8 | 3.925958e-8 | 1.26e-7 |
0.07 | 0.496585303791409 | 0.49658537342543 | 6.963402e-8 | 5.779747e-8 | 1.15 e-7 |
0.08 | 0.449328964117222 | 0.44932891673509 | 4.738213e-8 | 7.855516e-8 | 1.05 e-7 |
0.09 | 0.406569659740599 | 0.40657004499466 | 8.395048e-7 | 1.0158421e-7 | 5.70e-4 |
0.10 | 0.367879441171442 | 0.36788015032909 | 7.091575e-7 | 1.2697376e-7 | 5.17 e-4 |
Fig. 2. Plots of exact and computed solution for Problem 2
Problem 3. Consider a highly oscillatory test problem
Consider the second-order differential equation
\[
y” + \lambda^2 y = 0,\quad \lambda = 2,
\]
with initial conditions
\[
y(0)=1,\qquad y'(0)=2.
\]
Exact solution
The general solution is \(y(t)=A\cos(\lambda t)+B\sin(\lambda t)\). Using the initial conditions we get \(A=1\) and \(B=1\), hence
\[
y(t)=\cos(2t)+\sin(2t).
\]
(Opeyemi et al, 2024).
The new method is applied to solve problem 3 and then compared the result with
Opeyemietal., (2024) and Omole et al., (2018) as shown below:
Table 3: Comparison of absolute error for Problem 3( h = 1/100)
t | Exact solution | Computed solution | Error | Opeyemie tal., (2024), | Omole et al., (2018), |
0.01 | 1.019798673359911 | 1.019798673359911 | 0.0000e-0 | 2.2890e-13 | 3.409e-11 |
0.02 | 1.039189440847612 | 1.039189440847612 | 0.00000e-0 | 6.5075e-13 | 3.239e-11 |
0.03 | 1.058164546414649 | 1.058164546414649 | 0.00000e-0 | 6.2638e-13 | 3.465e-11 |
0.04 | 1.076716400271792 | 1.076716400129661 | 1.42131e-10 | 1.63787e-12 | 2.40e-13 |
0.05 | 1.094837581924854 | 1.094837586477904 | 4.55305e-9 | 3.88273e-12 | 1.780e-12 |
0.06 | 1.112520843142786 | 1.112520852389196 | 9.24641e-9 | 6.57699e-12 | 7.467e-11 |
0.07 | 1.129759110856874 | 1.129759124792946 | 1.393607e-8 | 9.70585e-12 | 3.904e-11 |
0.08 | 1.146545489989873 | 1.146545508672413 | 1.868254e-8 | 1.32541e-12 | 4.132e-11 |
0.09 | 1.162873266213946 | 1.162873287508854 | 2.129490e-8 | 1.72064e-11 | 1.197e-10 |
0.10 | 1.178735908636303 | 1.178735932535063 | 2.389876e-8 | 2.15469e-11 | 8.342e-11 |
Fig. 3. Plots of exact and computed solution for Problem 3
DISCUSSION OF RESULTS
The results obtained from the three test problems considered are summarized in Tables 1 – 3 and Figs. 1 – 3. In order to examine the usefulness and applicability of the new method, in table 1, we presented the numerical results and comparison of the absolute error with (Ogunware et al., (2023); Skwame et al., 2018) using the same order and values of h. It is very obvious that the new method is advantageous over other existing methods which solved the same problems in the literature. Table 2 shows the comparison of the results of our method with that of Opeyemi et al., (2024) and Akeremale et al., (2024) for problem 2. It is obvious that our method compares favorably with the existing methods despite their different methods. Table 3 shows the results generated for test problem 3 using the new method and that of Opeyemi et al., (2024) and Omole et al., (2018). The new method gives a minimal error. In general, figure 1-3 show that the implementation of the new scheme on three specific examples demonstrates its favorable comparison with exact solutions.
CONCLUSION
This paper has demonstrated how self-starting block methods for the solution of second order can be constructed by applying shift operator on two different linear multi-step formulae with low LTE. The continuous coefficients of the linear multi-step methods are derived using the undetermined coefficients method. Three problems in literature were used to show the efficiency of the methods and accuracy when compared to other methods. It was observed that the Block Methods for solving second order Ordinary Differential Equation using Shifted Linear Multi-step Formulae compare favorably well with some existing methods cited in the literature and Exact solution. The problem was solved using Matlab, not an in-build code but a code written by the authors.
REFERENCES
- Ajie, I. J., Durojaye, M. O., Utalor, I. K., & Onumanyi, P. A. (2020). Construction of a family of stable one-Block methods using Linear Multi-Step Quadruple.IOSRJournal of Mathematics, e- ISSN: 2278-5728, Vol (16), No. (4), pp. 01-13.
- Ajie, I. J., Utalor, I. K., & Onumanyi, P. A. (2019). A family of high order one-block methods for the solution of stiff initial value problems. Journal of Advances in Mathematics and Computer Science., 31(6), 1–14.
- Ajie, I. J., Utalor, I. K., &Durojaiye, M. O. (2019). Construction of stable high order one-block methods using multi-block triple. Journal of Advances in Mathematics and Computer Science., 32(5), 1–1
- Adesanya, A. O. (2011). Block methods for direct solution of general higher order initial
- value problems of ordinary differential equations, PhD Thesis, Federal University of
- Technology, Akure. (Unpublished).
- Adesanya, A. O., Anake, T. A., and Oghoyon, G. J. (2009). Continuous Implicit Method for the Solution of General Second Order Ordinary Differential Equations. Journal of Nigerian Association of Mathematical Physics, 15, 71–78.
- Adesanya, A. O., Anake, T. A. & Udoh, M. O. (2008). Improved continuous method
- for direct solution of general second order ordinary differential equations, Journal of Niger. Assoc. Math.Phys., 13, 59-62.
- AkeremaleC., Kuboye J.O.(2024). ‘’Hbride-Block Numerical Method for Ssolving Second Order Oordinary Differential Equations” International Journal of Computational Analysis. Vol (3), No. (2), pp. 26-35. https://www.researchgate.net/publication/377491717
- Awoyemi E.A., Adebile A.O., Adesanya A.O, Anake T.A, (2011) “Modified block method for the direct solution of second order ode” international Journal of Applied Mathematics and Computation. Vol 3(3) pp 181-188, 2011
- BadmusM,(2014) “ A new eight order implicit block algorithms for the direct solution of second order ordinary differential equations”, American Journal of Computational Mathematics, 376.https://doi.org/10.4236/ajcm.2014.44032.
- Badmus, A. M. & Yahaya, Y. A. (2009). An accurate uniform order 6 for the directsolution of general second order ordinary differential equations. The Pacific of Science and Technology.10 (2), 248-253
- Brugnano, L., and Trigiante, D. (1998). Solving Differential Problems by Multistep Initial and Boundary Value Methods. Amsterdam: Gordon and Breach Science Publishers.
- Fatunla, S. O. (1988). Numerical methods for initial value problems in ordinary differential equations, Academic press inc. Harcourt Brace Jovanovich Publishers, New York.
- Henrici, P.(1962). Discrete Variable Methods in ODE. New York: John Wiley & Sons.
- Lambert J. D.(1973), Computational methods in Ordinary Differential Equations, John Wiley and sons, New York.
- OgunwareG., Okafor F.M, Omole E.O. & Awoyemi F.C,(2023) “Direct Solution of second order ordinary differential equations with a one-step hybrid numerical model”, KIU Journal of Science, Engineering and Technology.Vol. 02 Issue 1, 45-52, ISSN: 1958-0641,https://doi.org/10.59568/KJSET-2023-2-1-07
- Okunnuga, S. O. (2008). One leg multi-step method for numerical integration of periodic second order initial value problems. N. A. M. P. 13, 63 – 68.
- Omar & Kuboye J.O.,(2015) “Application of order nine block method for solving second order ordinary differential equations directly”, Research Journal of Applied Sciences, Engineering and Technology.https://doi.org/10.19026/rjaset.11.1671.
- Omole E.O.& Ogunware B. G.(2018), “3- Point single hybrid block method (3PSHBM) for direct solution of general second order initial value problem of ordinary differential equations”, Asian Journal of Scientific Research and Reports 20 . https://org/10.9734/JSRR/2018/19862
- Onumanyi, Serisena U.W. & Jator S.N.,(1999) “Continuous finite difference approximation for solving differential equations”, International Journal of ComputerMathematics 72.https://doi.org/10.1080/00207169908804831.
- Opeyemi O. E, Catherine O. A” A fifth order block methods for solving second-order stiff ordinary differential equations using trigonometric functions and polynomial function as the basis function” Nigeria society of physical science.African Scientific Reports 3 (2024) 156
- Skwame Y, Sabo J, & Mathew M. (2019). The treatment of second order ordinary differential equations using equidistant one-step block hybrid. Asian Journal of Probability and Statistics. 5(3), 1-9.
- Utalor, K., Ajie, I. J., Durojaye, M.O, (2025). Construction of Block Hybrid Methods for the Solution of Singular Second Order Ordinary Differential Equation. Asian Research Journal of Mathematics,21(5),23–37. https://doi.org/10.9734/arjom/2025/v21i5923
- Vigo-Aguilar, J., & Ramos, H. (2006).Variable Stepsize Implementation of Multistep Methods for Y’(t)=f(x,y,y). Journal of Computational & Applied Mathematics,192,