J. Nig. Soc. Phys. Sci. 2 (2020) 160–165

Journal of the
Nigerian Society

of Physical
Sciences

Original Research

L-Stable Block Hybrid Numerical Algorithm for First-Order
Ordinary Differential Equations

B. I. Akinnukawea,∗, K. O. Mukab

aDepartment of Mathematics, University of Lagos, Lagos, Nigeria
bDepartment of Mathematics, University of Benin, Benin City, Nigeria

Abstract

In this work, a one-step L-stable Block Hybrid Multistep Method (BHMM) of order five was developed. The method is constructed for solving
first order Ordinary Differential Equations with given initial conditions. Interpolation and collocation techniques, with power series as a basis
function, are employed for the derivation of the continuous form of the hybrid methods. The discrete scheme and its second derivative are derived
by evaluating at the specific grid and off-grid points to form the main and additional methods respectively. Both hybrid methods generated are
composed in matrix form and implemented as a block method. The stability and convergence properties of BHMM are discussed and presented.
The numerical results of BHMM have proven its efficiency when compared to some existing methods.

Keywords: Second derivative, Stability, Hybrid, Block, Collocation techniques

Article History :
Received: 25 May 2020
Received in revised form: 13 July 2020
Accepted for publication: 14 July 2020
Published: 01 August 2020

c©2020 Journal of the Nigerian Society of Physical Sciences. All rights reserved.
Communicated by: J. O. Kuboye

1. Introduction

Most Mathematical models are formulated using ordinary
differential equations (ODEs) of different orders. However, math-
ematical models that are based on ODEs of order one occur in
several engineering, applied sciences and economics problems.
Some of the ODEs have been proven not to have closed-form
solution hence the need to develop numerical methods to pro-
vide approximate solutions to these problems. Consider the Ini-
tial Value Problem of the form

y′(x) = f (x, y(x)), y(x0) = y0 (1)

Several numerical methods have been developed by different
researchers to circumvent the Dahlquist’s order Barrier Theo-

∗Corresponding author tel. no:
Email address: iakinnukawe@unilag.edu.ng (B. I. Akinnukawe )

rem [1 - 3] that restricted the use of Linear Multistep Methods
of high order to solve (1). Some of such researchers are En-
right [4], Akinnukawe and Okunuga [5], Gear [6], Okunuga
[7], Cash [8], Akinfenwa et al. [9, 10], Adesanya et al. [11],
Ijezie and Muka [12] to mention but a few. High-order A-stable
and L-stable numerical methods are developed by incorporat-
ing off-step points, additional stages and/or employing higher
differentiation of the solution (Lambert [13]). Generating ap-
proximate solutions for (1) in instances where they exist and
are unique but are insoluble via analytical means is very impor-
tant because they help in the analysis and validation of models
in which they evolve from. One of the objectives of developing
numerical schemes for solving (1) is to obtain methods with
wider stability regions, better convergent rates and computa-
tional efficiency.

160



Akinnukawe & Muka / J. Nig. Soc. Phys. Sci. 2 (2020) 160–165 161

In this article, Continuous Hybrid Methods (CHM) of order 5
are derived through Interpolation and Collocation techniques
(Onumanyi et al. [14]). Both methods incorporate off-step
points and the second derivative of the scheme to produce the
discrete hybrid methods which are composed and implemented
as a block method (BHMM) to simultaneously produce approx-
imation at nodal points. The BHMM has the advantage of being
self-starting, L-stable in nature and possesses high accuracy be-
cause it implemented as a block method (see [15], [16], [17]).

2. Derivation of BHMM

This section describes the derivation of the k-step hybrid
multistep method of the form:

k∑
j=0

α jyn+ j + αvyn+v = h
k∑

j=0

β j fn+ j + hβv fn+v + h
2φkgn+k (2)

where h, n and v = 12 are the step size, grid index and off-
step point respectively while α j, β j, j = 0, 1, ..., k and φk are
parameters to be determined uniquely. An approximate solution
to (1) by the interpolating function

y(x) =
4k+1∑
j=0

b jt
j (3)

where b j, j = 0, 1, ..., 4k + 1 are the unknown coefficients. Im-
posing condition for the construction of the proposed class of
methods are:

y(xn+r ) = yn+r, r = 0, v (4)

y′(xn+r ) = fn+r, r = 0(v)k (5)

y′′(xn+r ) = gn+r, r = k (6)

Equations (4) - (6) will lead to a system of 4k + 2 equations.
These equations are solved simultaneously to obtain b j and the
values of b j are substituted into (3) to form the Continuous Hy-
brid Method (CHM) expressed in the form

y(t) =
k−v∑
j=0

α j(t)yn+ j + h
k∑

j=0

β j(t) fn+ j + h
2φk(t)gn+k (7)

where t = x−xn+k−1h and α j(t), β j(t), j = 0(v)k and φk(t) are the
continuous coefficients. The main scheme is generated evaluat-
ing CHM (7) at xn+k while the additional scheme is generated
from the second derivative of (7) at xn+v of the form

h2y′′(t) =
k−v∑
j=0

α j(t)yn+ j + h
k∑

j=0

β j(t) fn+ j + h
2φk(t)gn+k (8)

The developed main and additional schemes are combined
and implemented simultaneously as a block hybrid multistep
method BHMM for the numerical integration of IVPs (1).

Following the steps discussed above, the Continuous Hy-
brid Method (CHM) for k=1 is equation (7)

y(t) =
k−v∑
j=0

α j(t)yn+ j + h
k∑

j=0

β j(t) fn+ j + h
2φk(t)gn+k

where

 α0(t)
α 1

2
(t)

 = [ 1 0 −48023 128023 −120023 384230 0 48023 −128023 120023 −38423
]


t0

t1

t2

t3

t4

t5


(1)


β0(t)
β 1

2
(t)

β1(t)

 =


0 1 −13123
265
23 −

224
23

68
23

0 0 −12823
464
23 −

504
23

176
23

0 0 1923 −
89
23

128
23 −

52
23





t0

t1

t2

t3

t4

t5


(2)

[
φ1(t)

]
=
[

0 0 − 746
17
23 −

26
23

12
23

]


t0

t1

t2

t3

t4

t5


(3)

Interpolating (7) at x = xn+1 to generate the main method at
k = 1 becomes

yn+1 =
7yn
23

+
16yn+ 12

23
+

h fn
23

+
8h fn+ 12

23
+

6h fn+1
23

−
h2gn+1

46
(9a)

Equation (8), the second derivative of (7) for k = 1 is

h2y′′(t) =
k−v∑
j=0

α j(t)yn+ j + h
k∑

j=0

β j(t) fn+ j + h
2φk(t)gn+k

where

 α0(t)
α 1

2
(t)

 = [ −96023 768023 −1440023 768023960
23 −

7680
23

14400
23 −

7680
23

]
t0

t1

t2

t3

(4)

β0(t)
β 1

2
(t)

β1(t)

 =

−

262
23

1590
23 −

2688
23

1360
23

−
256
23

2784
23 −

6048
23

3520
23

38
23 −

534
23

1536
23 −

1040
23




t0

t1

t2

t3

 (5)
161



Akinnukawe & Muka / J. Nig. Soc. Phys. Sci. 2 (2020) 160–165 162

[
φ1(t)

]
=
[
−

7
23

102
23 −

312
23

240
23

]


t0

t1

t2

t3

 (6)
Interpolating(8) at x = xn+1/2, the additional method at k =

1 becomes

h2gn+ 12 =
240yn

23
−

240yn+ 12
23

+
31h fn

23
+

64h fn+ 12
23

+
25h fn+1

23
−

4h2gn+1
23

(9b)

The discrete hybrid methods (9a)and (9b) together forms
the One-step Block Hybrid Multistep Methods (BHMM). The
BHMM can be presented in a matrix block form as

A(1)Y$ = A(0)Y$−1 + hB(1) F$ + hB(0) F$−1 + h
2C(1)G$ (10)

where

Y$ =
 yn+ 12

yn+1

 ; Y$−1 =  yn− 12
yn

 ; F$ =  fn+ 12
fn+1

 ;
F$−1 =

 fn− 12
fn

 ; G$ =  gn+ 12
gn+1


The 2 by 2 matrices A(0), A(1), B(0), B(1), C(1), D(1) of the BHMM

(9) are defined as follows

A(1) =


240
23 0

−
16
23 1



A(0) =


0 24023

0 723



B(1) =


64
23

25
23

8
23

6
23



B(0) =


0 3123

0 123



C(1) =


−1 − 423

0 − 146



3. Analysis of BHMM

3.1. Order and Error Constant of the Method
Following Lambert [13] and Fatunla [18], a method was

proposed for finding the order p and error constant Wp+1 of the
block method (9) by first expanding y−, f− and g− functions by
Taylors series expansion about x and then comparing the coef-
ficients of h. It is established from our calculation that one-step
Block Hybrid Multistep Method have order and error constants
as p = (5, 5)T and Wp+1 = (

13
44160 ,

1
66240 )

T respectively where T
is transpose.

3.2. Zero Stability
A numerical method is said to be zero-stable if the roots

R j, j = 1, 2, . . . , N of the first characteristic polynomial ρ(R)
satisfies |R j| ≤ 1, j = 1, . . . , N and those roots with |R j| = 1 is
simple (see Lambert [3]). Applying the above conditions to the
derived block method, the first characteristic polynomial ρ(R) =
0 is calculated as

ρ(R) = det(RA(1) − A(0)) =
240
23

R(R − 1)

The BHMM is found to be zero-stable since ρ(R) = 0 satisfies
|R j| ≤ 1, j = 1, 2.

3.3. Convergence
According to Henrici [19], a numerical method converges

if it is consistent and zero-stable. Since BHMM (9) is of order
5 > 1, then it is consistent and we have established earlier that
the method satisfies the conditions of zero-stability. Therefore,
the block method (9) converges.

3.4. Stability of BHMM
Applying the BHMM to the test equation

y′ = λy, λ ≤ 0

we obtain

Y$ = Q(z)Y$−1, z = λh

where Q(z) is the amplification matrix given by

Q(z) =
A(0) + zB(0)

A(1) + zB(1) + z2C(1)
Q(z) has eigenvalues (ζ1,ζ2) = (0,ζ2). The dominant eigen-

value ζ2 is the stability function with real coefficient as

ζ2 =
10.4348 + 4.17391z + 0.652174z2 + 0.0434783z3

10.4348 − 6.26087z + 1.69565z2 − 0.26087z3 + 0.0217391z4

The stability function is used to plot the Region of Absolute
Stability (RAS) of the BHMM (see Figure 1). The proposed
method BHMM is L-stable since the RAS covers the entire left
plane of the graph (A-stable) and the limit of the stability func-
tion ζ2 is zero as z →∞

162



Akinnukawe & Muka / J. Nig. Soc. Phys. Sci. 2 (2020) 160–165 163

Figure 1. Stability Region of Block Hybrid Multistep Method

4. Numerical Results

The following problems are considered to examine the ac-
curacy and computational efficiency of BHMM. The computa-
tions were carried out using MATHEMATICA 9.0 software.
The following acronyms are used:

• ES = Exact Solution

• ESDMM = Error in SDMM

• EMHIRK = Error in MHIRK

• EBHMM = Error in BHMM

• EESDMM = Error in ESDMM

Problem 1: Consider the non-linear IVP

y′1 = λy1 + y
2
2

y′2 = −y2

with initial conditions

y1(0) =
−1
λ + 2

, y2(0) = 1,

where λ = 104 and the exact solution of the IVP is given as

y1(x) = −
e−2x

λ + 2
, y2(x) = e

−x

In Table 1, the numerical results obtained using BHMM
with h = 10−1 compares favourably with MHIRK method [10]
with h = 10−1 and is superior to that of Hojjati et al. [20] that
used h = 10−4.

Problem 2: Consider the following stiff problem arising
from Chemical Kinetic reactions in a Chemistry experiment.

y′1 = −0.013y1 − 1000y1y2 − 2500y1y3

y′2 = −0.013y1 − 1000y1y2

y′3 = −2500y1y3

with initial conditions

y1(0) = 0, y2(0) = 1, y3(0) = 1

The computed result at x = 2.0 from the new method BHMM
is compared with those of Ismail and Ibrahim (ESDMM [21]),
Hojjati et al. SDMM [20], Akinfenwa et al. MHIRK [10].
The step length h = 0.0125 was used for MHIRK method
and the new method BHMM and compared with Ismail and
Ibrahim ESDMM [21], Hojjati et al. SDMM [20] with step
length h = 0.001. The result obtained with the new method
BHMM is superior to others. See Table 2 for the numerical re-
sults.

Problem 3: Consider the non-linear system

y′1 = −2y1 + y2 + 2sin x

y′2 = 998y1 − 999y2 + 999(cos x − sin x)

with initial conditions

y1(0) = 2, y2(0) = 3, 0 ≤ x ≤ 10

with analytical solution given as

y1(x) = 2e
−x

+ sin x, y2(x) = 2e
−x

+ cos x

163



Akinnukawe & Muka / J. Nig. Soc. Phys. Sci. 2 (2020) 160–165 164

Table 1. Comparison of Methods for Problem 1

x yi ES ESDMM[20] EMHIRK[10] EBHMM
3 y1 −0.2478257E − 06 2.478147E − 11 3.06450E − 15 5.00564E − 16

y2 0.4978707E − 01 2.471093E − 06 3.07825E − 10 5.02813E − 11
5 y1 −0.4539085E − 08 3.450217E − 14 9.35475E − 17 1.52787E − 17

y2 0.6737946E − 02 2.304573E − 08 6.94326E − 11 1.13414E − 11
10 y1 0.3059023E − 12 3.456372E − 18 8.49412E − 21 3.75372E − 20

y2 0.3059023E − 04 3.150734E − 10 9.35666E − 13 1.52836E − 13

Table 2. Comparison of Methods for Problem 2

x yi ES EESDMM[21] ESDMM[20] EMHIRK[10] EBHMM
2.0 y1 −0.361693316929E − 5 0.82E − 10 0.52E − 13 0.110E − 13 2.919E − 15

y2 0.9815029948230 0.61E − 05 0.19E − 08 0.220E − 08 5.586E − 10
y3 1.018493388244 0.57E − 05 0.63E − 08 0.220E − 08 5.584E − 10

The numerical results of Problem 3 are shown in Table 3
using step size h = 10−3. The derived method integrated the
problem efficiently that the numerical solution is close to that
of the analytical solution.

Table 3. The absolute error for Problem 3
x yi ES EBHMM

0.25 y1 1.805005525 4.50751E − 14
y2 2.526513988 4.84057E − 14

0.5 y1 1.692486858 9.85878E − 14
y2 2.090643881 9.81437E − 14

1.0 y1 1.577229867 9.45910E − 14
y2 1.276061188 9.54792E − 14

2.0 y1 1.179967993 1.68310E − 13
y2 −0.145476270 1.68365E − 13

4.0 y1 −0.720171217 2.21378E − 13
y2 −0.617012343 2.23044E − 13

6.0 y1 −0.274457993 1.01363E − 13
y2 0.9651277911 1.01474E − 13

8.0 y1 0.9900291719 1.93401E − 13
y2 −0.1448291085 1.94650E − 13

10.0 y1 −0.5439303110 6.10623E − 13
y2 −0.8389807292 6.09068E − 13

Problem 4: A stiff system of Initial Value Problems

y′1 = −8y1 + 7y2

y′2 = 42y1 − 43y2

with initial conditions

y1(0) = 1, y2(0) = 8, x ∈ [0, 15]

with exact solution given as

y1(x) = 2e
−x
− e−50x, y2(x) = 2e

−x
+ 6 − e−50x

Problem 4 was integrated using the step size of h = 10−4 to
aid in comparing with other methods in literature as shown in
Tables 4 and 5. It is discovered that BHMM has better accuracy
than the others compared with.

Table 4. Numerical Results for Problem 4
x yi ES EBHMM
3 y1 9.9574136 × 10−2 2.68577 × 10−13

y2 9.9574136 × 10−2 2.65843 × 10−13

6 y1 4.9575044 × 10−3 1.68580 × 10−14

y2 4.9575044 × 10−3 1.80611 × 10−14

9 y1 2.4681961 × 10−4 7.57646 × 10−15

y2 2.4681961 × 10−4 5.43191 × 10−15

12 y1 1.2288424 × 10−5 2.10193 × 10−15

y2 1.2288424 × 10−5 2.54783 × 10−15

15 y1 6.1180464 × 10−7 2.29273 × 10−14

y2 6.1180464 × 10−7 1.87085 × 10−14

5. Conclusion

A one-step L-stable Block Hybrid Multistep Method (BHMM)
of order five was developed via interpolation and collocation
techniques. The BHMM has the advantage of being self-starting
and its L-stable in nature as displayed in figure 1. The block
method possesses high accuracy as shown in Tables (1) − (5)
where it was compared with some existing methods. The method
satisfies the zero-stability, consistency and convergence condi-
tions. BHMM has proved efficient for solving first-order Initial
Value Problems.

164



Akinnukawe & Muka / J. Nig. Soc. Phys. Sci. 2 (2020) 160–165 165

Table 5. A Comparison of absolute errors of methods for Problem 4

x Methods Error |y(tn) − yn|
5 BHMM 3.1999 × 10−14

TDBDF [22] 1.5472 × 10−2

TDMM [23] 1.5476 × 10−2

10 BHMM 4.642 × 10−15

TDBDF [22] 9.0808 × 10−5

TDMM [23] 9.0808 × 10−5

15 BHMM 1.8709 × 10−14

TDBDF [22] 6.1186 × 10−7

TDMM [23] 6.1186 × 10−7

Acknowledgments

The authors wish to thank the referees and editor for the
comments and suggestions to make this paper a success.

References

[1] G. G. Dahlquist, “A special stability problem for linear multistep meth-
ods”, BIT 3 (1963) 27.

[2] J. C. Butcher, The numerical methods for ordinary differential equations,
John Wiley and Sons Ltd, Chichester, 2008.

[3] J. D. Lambert, Computational methods for ordinary differential systems:
The Initial Value Problems, Wiley, Chichester, 1973.

[4] W. H. Enright, “Second derivative multistep methods for stiff ODEs”,
SIAM Journal of Numerical Analysis 11 (1974) 321.

[5] B. I. Akinnukawe & S. A. Okunuga, “A Seventh-order block integrator for
solving stiff systems”, Nigerian Journal of Mathematics and Applications
24 (2015) 67.

[6] C. W. Gear, “Hybrid methods for initial value problems in ordinary dif-
ferential equations”, SIAM Journal of Numerical Analysis 2 (1965) 69.

[7] S. A. Okunuga, “A fourth order composite two step method for stiff prob-
lems”, International Journal of Computer Mathematics 72 (1999) 39.

[8] J. R. Cash, “On the integration of stiff systems of ODEs using Extended
Backward Differentiation Formulae”, Numerische Mathematik 34 (1980)
235.

[9] O. A. Akinfenwa, B. Akinnukawe & S. B. Mudasiru, “A family continu-
ous third derivative block methods for solving stiff systems of first order
ODEs”, Journal of the Nigerian Mathematical Society 34 (2015) 160.

[10] O. A. Akinfenwa, S. A. Okunuga, B. I. Akinnukawe, U. P. Rufai & R.
I. Abdulganiy, “Multi-derivative hybrid implicit Runge-Kutta method for
solving stiff system of a first order differential equation”, Far East Journal
of Mathematical Sciences (FJMS) 106 (2018) 543.

[11] A. O. Adesanya, R. O. Onsachi & M. R. Odekunle, “New algorithm for
first order stiff initial value problem”, Fasciculi Mathematici 58 (2017) 2.

[12] A. P. Ijezie & K. O. Muka, “Modified SDBDF based on a non-zero root
of the second characteristics polynomial”, Journal of Science and Tech-
nology Research 1 (2019) 128.

[13] J. D. Lambert, Numerical methods for ordinary differential systems, John
Wiley, New York, 1991.

[14] P. Onumanyi, U. W. Sirisena & S. N. Jator. Continuous finite difference
approximations for solving differential equations, “International Journal
of Computational Mathematics”, 72 (1999) 15.

[15] W. E. Milne, Numerical solution of differential equations, John Wiley,
New York, 1953.

[16] D. Sarafyan, Multistep methods for the numerical solution of ODEs made
self-starting, Technical Report No. 495, Mathematics Research Center,
Madison, Wisconsin, 1965.

[17] L. F. Shampine & H. A. Watts, “Block implicit one-step methods”, Math-
ematics of Computation 23 (1969) 731.

[18] S. O. Fatunla, “Block methods for second order IVPs”, Internal Journal
of Computational Mathematics 41 (1991) 55.

[19] P. Henrici, Discrete Variable Methods in ODEs, John Wiley, New York,
1962.

[20] G. Hojjati, M. Y. Rahimi-Ardabili & S. M. Hosseini, “New second deriva-
tive multistep methods for stiff systems”, Applied Mathematical Mod-
elling 30 (2006) 466.

[21] G. A. Ismail & I. H. Ibrahim, “New efficient second derivative multistep
methods for stiff systems”, Applied Mathematical Modelling 23 (1978)
279.

[22] R. I. Okuonghae, M. N. O. Ikhile & J. Osemeke, “An off-step-point meth-
ods in multistep integration of stiff ODEs”, NMC Journal of Mathematical
Sciences 3 (2014) 731.

[23] A. K. Ezzeddine & G. Hojjati, “Third derivative multistep methods for
stiff systems”, International Journal of Nonlinear Science 14 (2012) 443.

165