J. Nig. Soc. Phys. Sci. 2 (2020) 134–140

Journal of the
Nigerian Society

of Physical
Sciences

Original Research

Linear Stability Analysis of Runge-Kutta Methods for Singular
Lane-Emden Equations

M. O. Ogunnirana,∗, O. A. Tayoa, Y. Harunab, A. F. Adebisia

aDepartment of Mathematical Sciences, Osun State University, Osogbo
bDepartment of Mathematics, Saadatu Rimi College of Education Kumbotso, Kano State

Abstract

Runge-Kutta methods are efficient methods of computations in differential equations, the classical Runge-Kutta method of order 4 happens to
be the most popular of these methods, and most times it is attached to the mind when Runge-Kutta methods are mentioned. However, there are
numerous forms of them existing in lower and higher orders of the classical method. This work investigates the linear stabilities and abilities of
some selected explicit members of these Runge-Kutta methods in integrating the singular Lane-Emden differential equations. The results obtained
established the ability of the classical Runge-Kutta method and why is mostly used in computations.

Keywords: Classical, Stability, Explicit, Lane-Emden equations, Ordinary Differential Equations

Article History :
Received: 13 April 2020
Received in revised form: 30 May 2020
Accepted for publication: 01 June 2020
Published: 01 August 2020

c©2020 Journal of the Nigerian Society of Physical Sciences. All rights reserved.
Communicated by: T. Latunde

1. Introduction

The most popular of the methods of Runge-Kutta is the clas-
sical Runge-Kutta method of order 4, ’classical’ as related to
the method obtained in the pre-computer era. However, many
other forms of Runge-Kutta methods have been derived by nu-
merous authors in the past. Existing methods of Runge-Kutta
are of orders 2, 3, and orders greater than 4. Methods of orders
greater than 4 are regarded as higher-order methods of Runge-
Kutta. The numerical solutions of singular Initial Value Prob-
lems (IVP) and Boundary Value Problems (BVP) of second-
order ordinary differential equations (ODEs) have been studied
in this work. Existing methods for the singular problems are

∗Corresponding author tel. no: +2347039104606
Email address: muideen.ogunniran@uniosun.edu.ng (M. O.

Ogunniran )

series, analytical methods and of recent non-series numerical
methods as obtained by various authors.
The general form of second-order singular differential equa-
tions is denoted as:

y′′(x) +
P(x)
Q(x)

y′(x) + f (x, y) = g(x, y); a ≤ x ≤ b. (1)

To solve equation (1), any of the conditions stated below need
to be imposed.

y(a) = α, y′(a) = β (2)

y(a) = α1, y
′(b) = β1 (3)

Equation (1) together with (2) are called Initial Value Problems
(IVPs) while equations (1) and (3) are called Boundary Value
Problems (BVPs).
Equation (1) is singular at Q(x) = 0, and f (x, y), g(x, y) are non-
linear continuous functions. It is well known that some of these

134



Ogunniran et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 134–140 135

problems have proved to be either difficult to solve or cannot
be solved analytically due to the singularity as the approximate
solutions lost their accuracy in the neighbourhood of the singu-
lar points. This discontinuous property has raised more interest
in many applied Mathematicians and Physicists in the studies
of this problem. Huta [1] developed the 8-stage Runge-Kutta
method of Order 6 in solving differential equation problems,
Qu and Agwarwal [2] employed collocation method for solv-
ing a class of singular non-linear two-point BVP. Koch, Peter
& Ewa [3] evaluated the approximate solution of the singular
IVP by implicit Euler method and finally used an acceleration
technique known as the iterated defect correction to improve
the approximations. Wazwaz [4, 5] presented series and exact
solution of Lane-Emden and Emden-Fowler type of problems
based on Adomian decomposition and modified Adomian de-
composition methods. A Simplified Derivation and Analysis
of Fourth Order Runge Kutta Method was presented by Musa,
Ibrahim, & Waziri [6]. Roul [7] presented a new efficient recur-
sive technique for solving singular boundary value problems
arising in various physical models. Abbas Al-Shimmary [8]
used the Runge-Kutta Method of order 6 to solve the initial
value problem of ordinary differential equations. Roul [9] also
presented a fourth-order B-spline collocation method and its er-
ror analysis for Bratu-type and Lane-Emden problems. A fast-
converging iterative scheme for solving a system of Lane–Emden
equations arising in catalytic diffusion reactions was presented
by Madduri and Roul [10]. Roul [11] presented a fast con-
verging iterative approach for the solution of doubly singular
BVP with derivative dependent source foundation. Ogunniran,
Haruna and Adeniyi [12] developed an efficient k−derivative
method for Lane-Emden equations and related stiff problems,
the paper focuses on the exploration of the possibility of a class
of multi-derivative methods for approximating its solution. A
new basic rational approximation method was developed for
solving singular initial value problems of ordinary differential
equations by Ogunniran & Adeniyi [13]. Ogunniran [14] also
developed a class of block multi-derivative numerical integra-
tor for singular advection equations. The solution of Emden-
Fowler equations was solved using the variational iteration method
by Olayiwola [15]. Olayiwola and Adegoke [16] presented
the approach of homotopy perturbation with Laplace transform
method in solving singular initial value problems. An opti-
mal sixth-order quartic B-spline collocation method was devel-
oped by Roul, Thula and Goura [17] for solving Bratu-type
and Lane-Emden–type problems. Roul, Madduri and Agar-
wal [18] developed a fast-converging recursive approach for
Lane–Emden type initial value problems arising in astrophysics.
Singh, Garg, Kanwar & Ramos [19] developed an Efficient Op-
timized Adaptive step-size Hybrid Block Method for Integrat-
ing Differential System. This study aims at the exploration of
the possibility of the classes of explicit Runge-Kutta methods

in the solution of singular Lane-Emden problems of ordinary
differential equations.

2. Method

The general m-stage Runge-kutta method is defined by

yr+1 − yr = hθ(xr, yr, h), (4)

θ(x, y, h) =
m∑

i=1

qiµi (5)

µ1 = f (x, y) (6)

µr = f (x + hpi, y + h
i−1∑
s=1

bisµs); i = 2, 3, · · · , m (7)

pi =
i−1∑
s=1

bis, i = 2, 3, · · · , m; h = xi − xi−1 (8)

It can be observed that an m-stage Runge-kutta method involves
m-function evaluations per step. Each of the functions µi(x, y, h),
i = 1, 2, 3, · · · , m may be interpreted as an approximation of the
derivative yi(x) and the function θ(x, y, h) as a weighted mean
of these approximations.
There is a great deal of tedious manipulation involved in de-
riving Runge-Kutta methods of the higher order. However, the
forms of Runge-Kutta methods in the scope of this work shall
be given in the following section.

(i) 2-stage Runge-kutta method of Order 2 (RK2)

yi+1 − yi = hµ2
µ1 = f (xi, yi)

µ2 = f (xi +
1
2 h, yi +

1
2 hµ1)

(9)

(ii) 3-stage Runge-kutta method of Order 3 (RK3)

yi+1 − yi =
1
6 h(µ1 + 4µ2 + µ3)

µ1 = f (xi, yi)
µ2 = f (xi +

1
2 h, yi +

1
2 hµ1)

µ3 = f (xi + h, yi − h(µ1 − 2µ2))

(10)

(iii) 4-stage Runge-kutta method of Order 4 (RK4)

yi+1 − yi =
1
6 h(µ1 + 2µ2 + 2µ3 + µ4)
µ1 = f (xi, yi)

µ2 = f (xi +
1
2 h, yi +

1
2 hµ1)

µ3 = f (xi +
1
2 h, yi +

1
2 hµ2)

µ4 = f (xi + h, yi + hµ3)

(11)

135



Ogunniran et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 134–140 136

(iv) 6-stage kutta-Nyström method of Order 5 (RK Nyström)

yi+1 − yi =
1

192 h(23µ1 + 125µ3 − 81µ5 + 125µ6)
µ1 = f (xi, yi)

µ2 = f (xi +
1
3 h, yi +

1
2 hµ1)

µ3 = f (xi +
2
5 h, yi +

1
25 h(4µ1 + 6µ2))

µ4 = f (xi + h, yi +
1
4 h(µ1 − 12µ2 + 15µ3))

µ5 = f (xi +
2
3 h, yi

+ 181 h(6µ1 + 90µ2 − 50µ3 + 8µ4))
µ6 = f (xi +

4
5 h, yi

+ 175 h(6µ1 + 36µ2 + 10µ3 + 8µ4))

(12)

(v) 8-stage Huta method of Order 6 (RK 8)

yi+1 − yi =
1

840 h(4µ1 + 216µ3 + 27µ4 + 272µ5 + 27µ6 + 216µ7 + 41µ8 )
µ1 = f (xi, yi )

µ2 = f (xi +
1
9 h, yi +

1
9 hµ1 )

µ3 = f (xi +
1
6 h, yi +

1
24 h(µ1 + 3µ2 ))

µ4 = f (xi +
1
3 h, yi +

1
6 h(µ1 − 3µ2 + 4µ3 ))

µ5 = f (xi +
1
2 h, yi +

1
8 h(−5µ1 + 27µ2 − 24µ3 + 6µ4 ))

µ6 = f (xi +
2
3 h, yi +

1
9 h(221µ1 − 9813µ2 + 867µ3 − 102µ4 + µ5 ))

µ7 = f (xi +
5
6 h, yi +

1
48 h(−183µ1 + 678µ2 − 472µ3 − 66µ4 + 80µ5 + 3µ6 ))

µ8 = f (xi + h, yi +
1
82 h(716µ1 − 2079µ2 + 1002µ3 + 834µ4 − 454µ5 − 9µ6 + 72µ7 ))

(13)

For purpose of completion and without loss of generality,
special Runge-kutta methods which were developed to
exhibit certain properties were also considered.

(vi) 5-stage Merson’s Method of Order 4 (RK Merson)

yi+1 − yi =
1
6 h(µ1 + 4µ4 + µ5)

µ1 = f (xi, yi)
µ2 = f (xi +

1
3 h, yi +

1
3 hµ1)

µ3 = f (xi +
1
3 h, yi +

1
6 h(µ1 + µ2))

µ4 = f (xi +
1
2 h, yi +

1
8 h(µ1 + 3µ3))

µ5 = f (xi + h, yi +
1
2 h(µ1 − 3µ3 + 4µ4))

(14)

(vii) 5-stage Scraton’s Method of Order 4 (RK Scraton)

yi+1 − yi = h
(

17
162 µ1 +

81
170 µ3 +

32
135 µ4 +

250
1377 µ5

)
µ1 = f (xi, yi)

µ2 = f (xi +
2
9 h, yi +

2
9 hµ1)

µ3 = f (xi +
1
3 h, yi +

1
12 h(µ1 + 3µ2))

µ4 = f (xi +
3
4 h, yi+

3
128 h(23µ1 − 81µ2 + 90µ3))

µ5 = f (xi +
9

10 h, yi+
9

10000 h(−345µ1 + 2025µ2 − 1224µ3 + 544µ4))

(15)

(viii) 6-stage England’s Method of Order 4 (RK England)

yi+1 − yi =
1
6 h(µ1 + 4µ3 + µ4)

µ1 = f (xi, yi)
µ2 = f (xi +

1
2 h, yi +

1
2 hµ1)

µ3 = f (xi +
1
2 h, yi +

1
4 h(µ1 + µ2))

µ4 = f (xi + h, yi − h(µ2 − 2µ3))
µ5 = f (xi +

2
3 h, yi+

1
27 h(7µ1 + 10µ2 + µ4))
µ6 = f (xi +

1
5 h, yi+

1
625 h(28µ1 − 125µ2 + 546µ3 + 54µ4 − 378µ5))

(16)

Figure 1. Stability Region for Runge-kutta Method of Order 2

2.1. Linear Stability

This is a behaviourial property related to h > 0. As in
most literature, the linear stability will be analyzed using the
Dalquist’s test

y′(t) = γy(t), <(γ) < 0. (17)

Applying methods (9)-(16) on (17), we have obtained a recur-
rence equation

yi = M(z)yi−1 (18)

where M(z = hγ) for each of the method is given in Table 1
below The contour for regions of absolute stability as obtained
from their characteristics equations are given in Figure 1 - 8.

3. Numerical Experiment

This section contains the numerical example considered and
their results presented in tables of absolute errors.

Example 1: Variable Coefficient Non-homogeneous Singular
IVP [3]

y′′(x) = −2x y
′(x) − n2cos(nx)

−
2
x nsin(nx); y(0) = 2, y

′(0) = 0

Theoretical Solution : y(x) = 1 + cos(nx)


(19)

* Absolute Error=|E xact solution − Numerical solution|, Nu-
merical solution is obtained using the Runge-kutta methods.

Example 2: Variable Coefficient Homogeneous Singular Initial
BVP [6]

(x2y′)′ = β(xy′ + y(α + β− 1))xα+β−2, 0 ≤ x ≤ 1.
y(0) = 1, y(1) = ex p(1)

Theoretical Solution = exp(x4)

(20)
136



Ogunniran et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 134–140 137

Table 1. Table of Characteristics Equations for Runge-kutta Methods

Method M(z), z = hλ
(9) 14 z

2 + 12 z + 1
(10) 16 z

3 + 12 z
2 + z + 1

(11) 124 z
4 + 16 z

3 + 12 z
2 + z + 1

(12) 1120 z
5 + 124 z

4 + 16 z
3 + 12 z

2 + z + 1
(13) z

8

483840 +
z7

4480 +
z6

720 −
179 z5
4032 −

3049 z4
3360 −

3119 z3
420 −

3119 z2
140 +

803 z
840 + 1

(14) 136 z
4 + 16 z

3 + 12 z
2 + z + 1

(15) 196 z
5 + 124 z

4 + 16 z
3 + 12 z

2 + z + 1
(16) 124 z

4 + 16 z
3 + 12 z

2 + z + 112

Table 2. Table of Absolute Errors for Example 1, n = 3 at different h

h x
y(x)

RK2 RK3 RK4 RK Nyström RK Merson RK Scraton RK England RK8

10−1
0.2 9.1982 × 10−2 9.0000 × 100 5.0392 × 10−2 4.5167 × 10−2 1.3355 × 10−1 7.61289 × 100 5.0392 × 10−2 1.1722 × 100

0.5 7.7896 × 10−2 9.0000 × 100 2.0474 × 10−2 7.5533 × 10−2 1.6385 × 10−1 8.3700 × 100 2.04738 × 10−2 1.2308 × 100

0.7 8.1030 × 10−2 9.0000 × 100 1.4680 × 10−2 8.1382 × 10−2 1.6970 × 10−1 7.4117 × 100 1.4680 × 10−2 1.2624 × 100

10−2
0.04 7.0148 × 10−4 9.0000 × 100 2.5148 × 10−4 7.0046 × 10−4 1.5808 × 10−3 9.7960 × 10−1 2.5148 × 10−4 1.2152 × 10−2

0.4 5.7493 × 10−4 9.0000 × 100 2.5188 × 10−5 9.2689 × 10−4 1.8072 × 10−3 2.7924 × 100 2.5188 × 10−5 3.9599 × 10−2

0.9 7.5147 × 10−4 9.0000 × 10+1 1.1203 × 10−5 9.4089 × 10−4 1.8212 × 10−3 1.200 × 10−2 1.1203 × 10−5 9.5724 × 10−2

10−3
0.09 4.7178 × 10−6 9.0000 × 100 1.1179 × 10−7 9.4074 × 10−6 1.8210 × 10−5 8.2736 × 10−1 1.1179 × 10−7 1.6996 × 10−3

0.19 4.8424 × 10−6 9.0000 × 10+1 5.2955 × 10−8 9.4662 × 10−6 1.8269 × 10−5 1.5242 × 100 5.2955 × 10−8 7.0497 × 10−3

0.69 6.7846 × 10−6 9.0000 × 10+1 1.4584 × 10−8 9.5046 × 10−6 1.8307 × 10−5 1.2185 × 100 1.4584 × 10−8 6.5798 × 10−2

Cmpt Time (s) 0.0025 0.0022 0.0035 0.0059 0.0039 0.0077 0.0065 0.021

Table 3. Table of Absolute Errors for Example 2 Using α = 2, β = 4

h x
y(x)

RK2 RK3 RK4 RK Nyström RK Merson RK Scraton RK England RK8

1
10

0.2 1.5913 × 10−3 1.5954 × 10−3 1.5745 × 10−3 1.5610 × 10−3 1.5706 × 10−3 7.0036 × 10−4 1.5745 × 10−3 1.5705 × 10−3

0.5 5.8226 × 10−2 5.8697 × 10−2 5.7050 × 10−2 5.7009 × 10−2 5.7014 × 10−2 2.5489 × 10−2 5.7050 × 10−2 5.7175 × 10−2

0.8 3.8263 × 10−1 3.8559 × 10−1 3.7177 × 10−1 3.7160 × 10−1 3.7166 × 10−1 5.7962 × 10−1 3.7177 × 10−1 3.7544 × 10−1

1
16

0.125 2.4357 × 10−4 2.4382 × 10−4 2.4258 × 10−4 2.4231 × 10−4 2.4234 × 10−4 1.5826 × 10−4 2.4258 × 10−4 2.4233 × 10−4

0.5 5.7534 × 10−2 5.7757 × 10−2 5.7026 × 10−2 5.7019 × 10−2 5.7020 × 10−2 2.5563 × 10−2 5.7026 × 10−2 5.7237 × 10−2

0.75 2.6863 × 10−1 2.8731 × 10−1 2.8304 × 10−1 2.8302 × 10−1 2.8303 × 10−1 3.7761 × 10−1 2.8304 × 10−1 2.8600 × 10−1

1
32

0.25 3.8056 × 10−3 3.8090 × 10−3 3.7977 × 10−3 3.7976 × 10−3 3.7976 × 10−3 1.1239 × 10−3 3.7977 × 10−3 3.8001 × 10−3

0.5625 9.0304 × 10−2 9.0405 × 10−2 9.0077 × 10−2 9.0076 × 10−2 9.0077 × 10−2 5.8501 × 10−2 9.0077 × 10−2 9.0639 × 10−2

0.875 5.5592 × 10−1 5.5618 × 10−1 5.5376 × 10−1 5.5538 × 10−1 5.5376 × 10−1 1.0939 × 100 5.5376 × 10−1 5.6323 × 10−1

Cmpt Time (s) 0.0023 0.0026 0.0042 0.0055 0.0033 0.0061 0.0055 0.011

137



Ogunniran et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 134–140 138

Figure 2. Stability Region for Runge-kutta Method of Order 3

Figure 3. Stability Region for Runge-kutta Method of Order 4

Figure 4. Stability Region for Kutta-Nyström Method

Example 3: Variable Coefficient Non-linear Homogeneous Sin-
gular IVP [3]

y′′(x) = −2x y
′(x) − y5(x); x ∈ [0, 1] y(0) = 1, y′(0) = 0

Theoretical Solution : y(x) = 1√
1+ x

2
3

(21)

Figure 5. Stability Region for Runge-kutta Method of Order 6

Figure 6. Stability Region for Merson’s Method

Figure 7. Stability Region for Scraton’s Method

Example 4: Van der Pol Singular Problem [11]

y′′ = y
′(1−y′2 )−y

µ
;

y(0) = 2, y′(0) = −23 +
10
81 µ−

292
2187 µ

2

−
1814

19683 µ
3; µ = 10−1.

(22)

The methods considered in this work were used to approximate
the problem over the interval [0, 0.55139] for h = 10−3.

138



Ogunniran et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 134–140 139

Table 4. Table of Absolute Errors for Example 3 Using different h

h x
y(x)

RK2 RK3 RK4 RK Nyström RK Merson RK Scraton RK England RK8

1
7

0.2857 2.4565 × 10−2 9.7182 × 10+6 2.7524 × 10−2 3.4496 × 10−2 4.1169 × 10−2 3.5739 × 10−2 2.7525 × 10−2 1.1722 × 10−1

0.5714 2.8265 × 10−2 9.7182 × 10+6 3.1850 × 10−2 3.8832 × 10−2 4.54985 × 10−2 3.9917 × 10−1 3.1851 × 10−2 1.2281 × 10−1

0.8571 1.0144 × 10−2 9.7182 × 10+6 1.4411 × 10−2 2.1093 × 10−2 2.7780 × 10−2 4.0945 × 10−1 1.4112 × 10−2 1.0723 × 10−1

1
14

0.2857 2.9836 × 10−2 1.2148 × 10+6 3.0700 × 10−2 3.2486 × 10−2 3.4150 × 10−2 2.1129 × 10−1 3.0701 × 10−2 5.4061 × 10−2

0.4286 3.4167 × 10−2 1.2148 × 10+6 3.5075 × 10−2 3.0686 × 10−2 3.8524 × 10−2 2.3129 × 10−1 3.5075 × 10−2 5.9065 × 10−2

0.6429 2.8969 × 10−2 1.2148 × 10+6 2.9944 × 10−2 3.1729 × 10−2 3.3393 × 10−2 2.4973 × 10−1 2.9944 × 10−2 5.5316 × 10−2

1
21

0.4286 3.4837 × 10−2 3.5993 × 10+5 3.5246 × 10−2 3.6042 × 10−2 3.6782 × 10−2 1.7714 × 10−1 3.5246 × 10−2 4.6596 × 10−2

0.6667 2.8266 × 10−2 3.5999 × 10+5 2.8707 × 10−2 2.9504 × 10−2 3.0243 × 10−2 1.9637 × 10−1 2.8707 × 10−2 4.1652 × 10−2

0.9048 9.3712 × 10−3 3.5993 × 10+5 9.8446 × 10−3 1.0641 × 10−2 1.1138 × 10−2 1.9780 × 10−1 9.8446 × 10−3 2.4823 × 10−2

Cmpt Time (s) 0.0025 0.0021 0.0035 0.0059 0.0045 0.0077 0.0051 0.021

Figure 8. Stability Region for England’s Method

Table 5. Numerical Results for Example 4 at x = 0.55139

h Method y(x) Cmpt. Time (s)

RK2 0.625821430245524 0.0025
RK3 0.625821606669667 0.0019
RK4 0.625821602543592 0.0065

10−3 RK Nyström 0.625821602551168 0.018
RK Merson 0.625821602551009 0.0028
RK Scarton -0.031597285014076 0.025
RK England 0.625821602543570 0.0021

RK 8 0.642280546551478 0.063

4. Discussion of Results and Conclusion

It is worthy to note that numerical methods were programmed
via MATLAB 9.2 version on a personal computer with the fol-
lowing specifications:

• System name- HP Pavilion x360 Convertible

• Processor- Intel(R) Core(TM) i3-7100U CPU @ 2.40GHz

• Installed memory (RAM)- 8.00GB

• System Type- 64-bits Operating System, x64-based pro-
cessor

• Operating system- Windows 10

Cmpt time is the computer computation time measured in
seconds(s). This paper presents the performances of different
orders of Runge-Kutta methods in the integration of singular
Lane-Emden equations. From Tables 2-4, It could be observed
that the Runge-Kutta method of order 4 performs the best in
terms of accuracy. The methods of order 3, Scraton’s, and order
8 were found to have failed in the solution of singular problems
in ordinary differential equations. The England’s method per-
forms more satisfactorily as it shows some competitive strength
against the order 4 method. It could be concluded that the
Runge-Kutta method of order 4 outperforms all other methods
under consideration and this property makes it the most stable
and accurate of the Runge-Kutta methods.

Acknowledgments

The authors thank the referees and editor for their valuable
comments and suggestions.

References

[1] A. Huta, “Contribution à la formule de sixieme ordre dans la methode de
Runge-kutta-Nyström”, Acta Fac Rerum Natu Univ. Comenian Math. 2
(1957) 21.

[2] R. Qu & R. P. Agarwal, “A collocation method for solving a class of sin-
gular non-linear two-point boundary value problems”, Journal Computer
Application Mathematics 83 (1997) 147.

[3] O. Koch, K. Peter & B. W. Ewa, “The Implicit Euler method for the nu-
merical solution of singular initial value problems”, Applied Numerical
Mathematics 34 (2000) 231

[4] A. M. Wazwaz, “A new algorithm for solving differential equation of
Lane-Emden type”, Applied Mathematics and Computation 118 (2001)
287.

[5] A. M. Wazwaz, “Adomian decomposition method for a reliable treatment
of the Emden-Fowler equation”, Applied Mathematics and Computation
61 (2005) 543.

[6] H. Musa, S. Ibrahim & M. Y. Waziri, “A Simplified Derivation and Anal-
ysis of Fourth Order Runge Kutta Method”, International Journal of Com-
puter Applications 9 (2010) 51.

[7] P. Roul, “A new efficient recursive technique for solving singular bound-
ary value problems arising in various physical models”, The European
Physical Journal Plus 131 (2016), doi.org/10.1140/epjp/i2016-16105-8.

139



Ogunniran et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 134–140 140

[8] A. Abbas Al-Shimmary, “Solving Initial Value Problem Using Runge-
kutta 6th Order Method”, ARPN Journal of Engineering and Applied Sci-
ences 12 (2017) 3953.

[9] P. Roul & K. Thula, “A fourth order B-spline collocation method and its
error analysis for Bratu-type and Lane-Emden problems”, International
Journal of Computer Mathematics 96 (2017) 85.

[10] H. Madduri, & P. Roul, “A fast-converging iterative scheme for solving a
system of Lane–Emden equations arising in catalytic diffusion reactions”,
Journal of Mathematical Chemistry 57 (2018) 570.

[11] P. Roul, “Doubly singular boundary value problems with derivative de-
pendent source function: A fast converging iterative approach”, Journal
of Math Methods Application Science 42 (2018) 354.

[12] M. O. Ogunniran, Y. Haruna, & R. B. Adeniyi, “Efficient k-derivative
Methods for Lane-Emden Equations and Related Stiff Problems”, Nige-
rian Journal of Mathematics and Applications 28 (2019) 1.

[13] M. O. Ogunniran & R. B. Adeniyi, “New Basic Rational Approximation
Method for Solving Singular Initial Value Problems of Ordinary Differen-
tial Equations”, ABACUS (Mathematics Science Series) 46 (2019) 119.

[14] M. O. Ogunniran, “A Class of Block Multi-Derivative Numerical Integra-
tor for Singular Advection Equations”, Journal of the Nigerian Society of

Physical Sciences 1 (2019) 62.
[15] M. O. Olayiwola, “Solutions of Emden-Fowler Type Equations by Varia-

tional Iteration Method”, Cankaya University Journal of Science and En-
gineering 16 (2019) 1.

[16] M. O. Olayiwola & A. Adegoke, “Homotopy Perturbation with Laplace
Transforms Method for Solving Singular Initial Value Problems” Nige-
rian Journal of Mathematics and Applications 28 (2019) 65.

[17] P. Roul, K. Thula & V. Goura, “An optimal sixth-order quartic B-spline
collocation method for solving Bratu-type and Lane-Emden–type prob-
lems”, Journal of Computational and Applied Mathematics 359 (2019)
182.

[18] P. Roul, H. Madduri & R. Agarwal, “A fast-converging recursive approach
for Lane-Emden type initial value problems arising in astrophysics”, Jour-
nal of Math Methods Application Science 42 (2019) 2613.

[19] G. Singh, A. Garg, V. Kanwar & H. Ramos, “An Efficient Optimized
Adaptive Step-size Hybrid Block Method for Integrating Differential Sys-
tems”, Applied Mathematics and Computation, 362 (2019) 124567.

[20] J. D. Lambert, Computational Methods in Ordinary Differential Equa-
tions, John Wiley and Sons, New York, 1991.

140