J. Nig. Soc. Phys. Sci. 4 (2022) 265–280

Journal of the
Nigerian Society

of Physical
Sciences

Approximate Analytical Solution of Fractional Lane-Emden
Equation by Mittag-Leffler Function Method

Richard Olu Awonusika∗, Oluwaseun Akinlo Mogbojuri

Department of Mathematical Sciences, Adekunle Ajasin University, P.M.B. 001, Akungba Akoko, Ondo State, Nigeria

Abstract

The classical Lane-Emden differential equation, a nonlinear second-order differential equation, models the structure of an isothermal gas sphere
in equilibrium under its own gravitation. In this paper, the Mittag-Leffler function expansion method is used to solve a class of fractional Lane-
Emden differential equation. In the proposed differential equation, the polytropic term f (y(x)) = ym(x) (where m = 0, 1, 2, . . . is the polytropic
index; 0 < x ≤ 1) is replaced with a linear combination f (y(x)) = a0 + a1y(x) + a2y2(x) + · · · + amym(x) + · · · + aN yN (x), 0 ≤ m ≤ N, N ∈ N0.
Explicit solutions of the fractional equation, when f (y) are elementary functions are presented. In particular, we consider the special cases of
the trigonometric, hyperbolic and exponential functions. Several examples are given to illustrate the method. Comparison of the Mittag-Leffler
function method with other methods indicates that the method gives accurate and reliable approximate solutions of the fractional Lane-Emden
differential equation.

DOI:10.46481/jnsps.2022.687

Keywords: Fractional differential equation, Lane-Emden differential equation, Mittag-Leffler function, Elementary functions, Series solution,
Eigenfunction expansion method

Article History :
Received: 02 March 2022
Received in revised form: 07 April 2022
Accepted for publication: 10 April 2022
Published: 29 May 2022

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

1. Introduction

Fractional calculus is a generalization and extension of clas-
sical calculus to non-integer orders. In recent times, fractional
calculus has attracted the attention of researchers in several ar-
eas including mathematics, physics, biology, chemistry, engi-
neering, economics and psychology ([1], [2], [3], [4], [5], [6],
[7], [8], [9]). Several definitions of fractional calculus have
been formulated by researchers. The most common and widely
used definitions are the Caputo, Grünwald-Letnikov and

∗Corresponding author tel. no:
Email address: richard.awonusika@aaua.edu.ng (Richard Olu

Awonusika)

Riemann-Liouville derivatives ([5], [8], [10]). The Grünwald-
Letnikov derivative is mostly limited to numerical algorithms.
The Riemann-Liouville fractional derivative has certain limi-
tations and so becomes unsuitable in modeling some real-life
phenomena since it requires the definition of fractional order
initial conditions which have no physical meaningful explana-
tion yet ([8]). The main advantage of the Caputo derivative
is that it takes on the same form of initial conditions for the
integer-order differential equations.

Recently, the Lane-Emden equation has been investigated
by several researchers due to its significant applications in math-
ematical physics and astrophysics ([11]). The classical Lane-
Emden equation, first introduced in 1870 by Lane and studied

265



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 266

in 1907 by Emden is of the form:

y′′(x) +
2
x

y′(x) + f (y(x)) = 0 (1)

with initial conditions

y(0) = A, y′(0) = B. (2)

Equation (1) with various forms of f (y(x)) has been used to
model several phenomena such as the theory of stellar struc-
ture, thermal explosions, the thermal behavior of a spherical
cloud of gas, isothermal gas spheres and thermionic currents
([12], [13]). Various methods have been presented by several
authors to obtain the solution of the initial value problem (1)
– (2). The Adomian decomposition method was employed in
[11] to investigate the initial value problem (1) – (2). Solutions
of the Lane-Emden equation have also been obtained via the
series method (see e.g., [14], [15]). The series solutions ob-
tained in [15] were compared with the results obtained using
the homotopy perturbation method. A numerical algorithm was
developed by Vanani and Aminataei [16] for solving the Lane-
Emden equation. Ogunniran et al. in [17] investigated the linear
stabilities of some explicit members of Runge-Kutta methods in
integrating the Lane-Emden equation.

However, the standard Lane-Emden differential equation is
the one given precisely by the initial value problem ([14])

y′′(x) +
2
x

y′(x) + ym(x) = 0, m ≥ 0; x > 0

y(0) = 1, y′(0) = 0,
(3)

where m is the polytropic index and y = y(x) is the polytrope.
Clearly, the equation (3) is linear when m = 0, 1 and as a result
the analytical solutions of the corresponding equations are real-
isable in closed forms. By extension it is mentioned in [14] that
a closed form solution is also possible for m = 5. For numer-
ical solutions of second order ordinary differential equations,
see [18] and [19].

As a result of the significant importance of fractional calcu-
lus in modeling real-life phenomena accurately, the fractional
Lane-Emden equation has been formulated and studied by re-
searchers in very recent times. By using the collocation method,
a numerical solution of (3) was obtained in [20]. Some other
researchers have also sought numerical solutions to the frac-
tional Lane-Emden equation (see, e.g., [21], [22]). Approxi-
mate solutions based on orthonormal Bernoulli’s polynomials
method, Homotopy-Adomian decomposition method and the
series expansion method were treated in [23], [24] and [25] re-
spectively. Other treatise on the fractional Lane-Emden equa-
tion can be found in [26], [27], [28], [29] and [30]. In [28, Sub-
section 5.2.2], the authors considered the power series solutions
of the fractional Lane-Emden equation with the polytropic term
ym with the fractional derivative described in the Caputo sense,
while Malik and Mohammed [25] presented approximate so-
lutions of fractional Lane-Emden equation using conformable
Homotopy–Adomian decomposition method and conformable
residual power series method.

Arafa et al. in [31] used the Mittag-Leffler function method
to solve a simple fractional differential equation of the form

Dαy = Ay2, 0 < α ≤ 1, (4)

with the fractional derivative described in the Caputo sense. In
this paper, the Mittag-Leffler function expansion method is em-
ployed to solve a class of fractional Lane-Emden initial value
problem whose polytropic term f (y(x)) = ym(x) (where m =
0, 1, 2, . . . is the polytropic index; 0 < x ≤ 1) is replaced
with a finite series f (y(x)) = a0 + a1y(x) + a2y2(x) + · · · +
amym(x) + · · · + aN yN (x), 0 ≤ m ≤ N, N ∈ N0, with the frac-
tional derivative described in the Caputo sense. The analytical
solutions of the corresponding fractional equations, when f (y)
are elementary functions, are presented, from which several ex-
amples of fractional Lane-Emden equations are given. In par-
ticular, we consider the special cases where f (y) are trigono-
metric, hyperbolic and exponential functions. Comparison of
the Mittag-Leffler function method with other methods shows
that the method is reliably capable of solving analytically, non-
linear fractional Lane-Emden differential equations, and by ex-
tension, one can apply the method to solve several nonlinear
fractional Lane-Emden equations when the functions f (y) are
given by other special functions. The motivation for the forms
of the nonlinear function f (y) considered in this paper arises
from the need to address those situations where the expansion
coefficients a0, a1, . . . , aN−1 are not identically zero.

2. Caputo Fractional Derivative and Its Basic Properties

This section presents the definition and basic properties of
Caputo derivative needed. For further discussions on fractional
calculus and fractional differential equations, see, e.g., [5], [28],
[32], [33], [34], [35], [36], [37], [38], [39], [40].

Definition 2.1. For α > 0, the Caputo fractional derivative of
order α is defined as follows (m ≥ 1) :

c Dα f (x) =Jm−αDm f (x)

=


1

Γ(m−α)

∫ x
0

(x − ξ)m−α−1 f (m)(ξ) dξ, m − 1 < α < m,

dm
d xm f (x) α = m,

(5)

where

J
α f (x) :=

1
Γ(α)

∫ x
0

(x − ξ)α−1 f (ξ) dξ (6)

is the Riemann-Liouville fractional integral of order α. Here
f (m)(x) = d

m f (x)
d xm .

Remark 2.1. The Caputo derivative of the constant function C
is zero, i.e c DαC = 0.

266



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 267

Corollary 2.1. For m − 1 < α < m, we have

c Dαxβ =



0, β ∈ {0, 1, · · · , m − 1} ,

Γ(ν+1)
Γ(β−α+1) x

β−α β > m − 1

non existing otherwise.

(7)

The following properties of the Caputo derivative hold ( [8],
[10]).

Lemma 2.1. Suppose m − 1 < α < m; n, m ∈ N; α ∈ R; λ,σ ∈
C. Let f (x) and g(x) be such that c Dα f (x) and c Dαg(x) exist.
We have that

(a) lim
α→n

c Dα f (x) = f (n)(x)

(b) lim
α→n−1

c Dα f (x) = f (n−1)(x) − f (n−1)(0)

(c) c Dα(λ f (x) + σg(x)) = λ c Dα f (x) + σ c Dαg(x)

(d) c DαDm f (x) = c Dα+m f (x) , Dm c Dα f (x).

3. Mittag-Leffler Function

In this section we examine the definition and properties of
the Mittag-Leffler function as eigenfunctions of fractional deriva-
tives ([1]). The Mittag-Leffler function plays an important role
in the theory of fractional calculus. Just as the exponential func-
tion naturally comes out of the solutions of classical differential
equations, the Mittag-Leffler function plays an analogous role
in the solution of fractional differential equations. It is thus a
generalization of the exponential function.

Definition 3.1. The Mittag-Leffler function is defined by:

Eα(z) =
∞∑

k=0

zk

Γ(αk + 1)
, α > 0,α ∈ R, z ∈ C. (8)

Definition 3.2. The Mittag-Leffler function can also be repre-
sented in two arguments α and β as

Eα,β(z) =
∞∑

k=0

zk

Γ(αk + β)
, α,β > 0,α,β ∈ R, z ∈ C; (9)

with Eα,1 = Eα.

In particular, the significance of the Mittag-Leffler function
is in the fact that it serves as the eigenfunction of the Caputo
and Riemann-Liouville derivatives in fractional calculus ([1]).
It is easy to see that ([8], [41])

E1,1(z) = E1(z) = e
z, E1,2(z) =

ez − 1
z

, E2,1(z
2) = cosh z, z ∈ C.

E2,1(−z
2) = cos z, E2,2(z

2) =
sinh z

z
, E2,2(−z

2) =
sin z

z
, z ∈ C.

Corollary 3.1. For α,β ∈ R, z ∈ C and m ∈ N, the following
recurrence relations hold for the Mittag-Leffler function ([41]) :

(a) Eα,β(z) = zEα,α+β(z) +
1

Γ(β)

(b) Eα,β(z) = βEα,β+1(z) + αz
d
dz Eα,β(z) = βEα,β+1(z)

(c)
(

d
dz

)m
[zβ−1 Eα,β(zα)] = zβ−m−1 Eα,β−m(zα).

4. Fractional Lane-Emden Equations and Their Approxi-
mate Solutions

In this section, we apply the Mittag-Leffler function dis-
cussed in Section 3 to obtain approximate solutions to the frac-
tional Lane-Emden initial value problem

c Dαy(x) +
ω

xα−β
c Dβy(x) + f (y(x)) = 0

y(0) = A, y′(0) = 0, A ∈ [0, 1],
(10)

with 0 < x ≤ 1, 0 < β ≤ 1, 1 < α ≤ 2, ω ∈ R. Here
c Dµ denotes the Caputo fractional differential operator of order
µ, µ > 0. The resulting solution is given as the Mittag-Leffler
function. The function f (y(x)) will be chosen in such a way that
it can be expanded in a power series, in particular, the Maclaurin
series.

Towards this end, let the function f (y(x)) be given by the
power series

f (y(x)) =
∞∑

m=0

amy
m(x), (11)

where am, m = 0, 1, 2, 3, . . . , are the expansion coefficients
(constants). For computation reasons, one is interested in the
approximation

fN (y(x)) =
N∑

m=0

amy
m(x), aN , 0. (12)

We proceed by giving the following result that will be needed
in the sequel.

Proposition 4.1. For 1 < α ≤ 2, 0 < x ≤ 1, the function
fN (y(x)) admits the power series expansion

FαN (x) := fN (y(x)) =
∞∑
`=0

Cα`,N x
α`, N ∈ N0, (13)

where the constants Cα
`,N are given by

Cα`,N =
N∑

m=0

amB
α,N
`,m . (14)

Here the numbers Bα,N
`,m , m = 3, 4, 5, . . . are given by the finite

series

Bα,N
`,m =

∑̀
pm−1 =0

pm−1∑
pm−2 =0

· · ·

p2∑
p1 =0

Aα`−pm−1,N A
α
pm−1−pm−2,N

· · ·Aαp1,N, (15)

with the Mittag-Leffler expansion coefficients given by

Aα`,N =
b`N

Γ(α` + 1)
, b`N are constants. (16)

In particular, the special cases Bα,N
`,m , m = 0, 1, 2, are given re-

spectively by

Bα,N
`,0 =


1, ` = 0,

0, ` = 1, 2, · · · ,

Bα,N
`,1 = A

α
`,N, B

α,N
`,2 =

∑̀
p1 =0

bp1N b
`−p1
N

Γ(αp1 + 1)Γ(α` −αp1 + 1)
. (17)

267



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 268

Proof. By definition, we use the Mittag-Leffler function

y(x) = Eα(bN x
α) =

∞∑
k=0

bkN x
αk

Γ(αk + 1)
(18)

in (12) to see that

FαN (x) := fN (y(x)) =
N∑

m=0

am

 ∞∑
k=0

bkN x
αk

Γ(αk + 1)

m

=

N∑
m=0

amS
α
m,N (x), (19)

where we have written (m ≥ 1)

Sα0,N = 1, S
α
m,N (x) := y

m(x) =

 ∞∑
k=0

bkN x
αk

Γ(αk + 1)

m . (20)
Let Aαk,N be given by

Aαk,N =
bkN

Γ(αk + 1)
. (21)

Then clearly,

Sα1,N (x) =
∞∑
`=0

Aα`,N x
α`

=

∞∑
`=0

Bα,N
`,1 x

α`. (22)

Furthermore, the classical Cauchy product tells us that

Sα2,N (x) =

 ∞∑
k=0

Aαk,N x
αk

2 = ∞∑
m=0

Aαm,N x
αm

∞∑
n=0

Aαn,N x
αn

=

∞∑
`=0

Bα,N
`,2 x

α`, (23)

where

Bα,N
`,2 =

∑̀
p1 =0

Aαp1,N A
α
`−p1,N

. (24)

Also, in a similar way, we see that

Sα3,N (x) =

 ∞∑
k=0

Aαk,N x
αk

3

=

 ∞∑
k=0

Aαk,N x
αk

2 ∞∑
n=0

Aαn,N x
αn. (25)

That is, one has

Sα3,N (x) = S
α
2,N (x)S

α
1,N (x)

=

 ∞∑
m=0

Bα,Nm,2 x
αm


 ∞∑

m=0

Bα,Nn,1 x
αn

 = ∞∑
`=0

Bα,N
`,3 x

α`, (26)

where

Bα,N
`,3 =

∑̀
p2 =0

Bα,Np2,2B
α,N
`−p2,1

=
∑̀
p2 =0

p2∑
p1 =0

Aαp1,N A
α
p2−p1,N

Aα`−p2,N. (27)

Similarly, by taking more products, we see that

Sα4,N (x) =

 ∞∑
k=0

Aαk,N x
αk

4

=Sα3,N (x)S
α
1,N (x) =

 ∞∑
m=0

Bα,Nm,3 x
αm


 ∞∑

n=0

Bα,Nn,1 x
αn


=

∞∑
`=0

∑̀
p3 =0

Bα,Np3,3B
α,N
`−p3,1

xα` =
∞∑
`=0

Bα,N
`,4 x

α`, (28)

where

Bα,N
`,4 =

∑̀
p3 =0

p3∑
p2 =0

p2∑
p1 =0

Aαp1,N A
α
p2−p1,N

Aαp3−p2,N A
α
`−p3,N

. (29)

Continuing in this way, we obtain

Sαm,N (x) =S
α
m−1,N (x)S

α
1,N (x)

=

 ∞∑
k=0

Bα,Nk,m−1 x
αk


 ∞∑

n=0

Bα,Nn,1 x
αn


=

∞∑
`=0

∑̀
pm−1 =0

Bα,Npm−1,m−1B
α,N
`−pm−1,1

xα` =
∞∑
`=0

Bα,N
`,m x

α`, (30)

where (m = 3, 4, 5, . . . )

Bα,N
`,m =

∑̀
pm−1 =0

pm−1∑
pm−2 =0

· · ·

p2∑
p1 =0

Aα`−pm−1,N A
α
pm−1−pm−2,N

· · ·Aαp1,N. (31)

Upon inserting (30) – (31) into (19), one gets

FαN (x) =
N∑

m=0

am
∞∑
`=0

Bα,N
`,m x

α`
=

∞∑
`=0

Cα`,N x
α`, (32)

where

Cα`,N =
N∑

m=0

amB
α,N
`,m (33)

and we obtain the result as required.

With the Proposition 4.1 in place we now present the main
result of this paper.

Theorem 4.1. For 0 < x ≤ 1, 0 < β ≤ 1, 1 < α ≤ 2, ω ∈
R, N ∈ N0; the fractional initial value problem

c Dαy(x) +
ω

xα−β
c Dβy(x) +

N∑
m=0

amy
m(x) = 0

y(0) = A, y′(0) = 0, A ∈ [0, 1],

(34)

268



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 269

admits a series solution given by the Mittag-Leffler function

y(x) =
∞∑
`=0

B`N (α,β; ω)
Γ(α` + 1)

xα`, (35)

where

B
0
N (α,β; ω) := b

0
N = A, (36)

B
`+1
N (α,β; ω) := b

`+1
N

= −
Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα

`,N

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
, ` ≥ 0. (37)

The coefficients Cα
`,N, ` ≥ 0, are as given in Proposition 4.1.

Proof. The starting point is to substitute the Mittag-Leffler func-
tion (18) in the equation (34) to evaluate the Caputo fractional
derivatives c Dαy(x) and xβ−α c Dβy(x). Towards this end, using
Corollary 2.1, we have that

c Dαy(x) =
∞∑

k=0

bkN
Γ(αk + 1)

c Dαxαk =
∞∑

k=1

bkN
Γ(αk + 1 −α)

xαk−α

=

∞∑
k=0

bk+1N
Γ(αk + 1)

xαk (38)

and similarly one sees that

xβ−α c Dβy(x) =xβ−α
∞∑

k=0

bkN
Γ(αk + 1)

c Dβxαk

=

∞∑
k=0

bk+1N
Γ(αk + α + 1 −β)

xαk. (39)

Substituting the series (38), (39) in (34) and then applying Propo-
sition 4.1 gives

∞∑
`=0

b`+1N
Γ(α` + 1)

xα` +
∞∑
`=0

ωb`+1N
Γ(α(` + 1) −β + 1)

xα`

+

∞∑
`=0

Cα`,N x
α`

= 0, (40)

which implies that

∞∑
`=0

 b`+1N
Γ(α` + 1)

+
ωb`+1N

Γ(α(` + 1) −β + 1)
+ Cα`,N

 xα` = 0. (41)
Comparing the coefficients of xα`, ` = 0, 1, 2, · · · on both sides
of (41) gives

b`+1N

(
1

Γ(α` + 1)
+

ω

Γ(α(` + 1) −β + 1)

)
= −Cα`,N

and as a result,

b`+1N = −
Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα

`,N

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
, (42)

where the coefficients Cα
`,N are as given in Proposition 4.1.

Remark 4.1. For computation purposes, in addition to the first
two coefficients Bα,N

`,m , m = 1, 2, given in Propostion 4.1, we
present explicitly the coefficients Bα,N

`,m , 3 ≤ m ≤ 20:

Bα,N
`,3 =

∑̀
p2 =0

p2∑
p1 =0

bp1N b
p2−p1
N b

`−p2
N

Γ(αp1 + 1)Γ(αp2 −αp1 + 1)Γ(α` −αp2 + 1)

Bα,N
`,4 =

∑̀
p3 =0

p3∑
p2 =0

p2∑
p1 =0

bp1N b
p2−p1
N b

p3−p2
N b

`−p3
N

Γ(αp1 + 1)Γ(αp2 −αp1 + 1)Γ(αp3 −αp2 + 1)Γ(α` −αp3 + 1)
...

Bα,N
`,10 =

∑̀
p9 =0

p9∑
p8 =0

· · ·

p2∑
p1 =0

b`−p9N b
p9−p8
N · · ·b

p1
N

Γ(α` −αp9 + 1)Γ(αp9 −αp8 + 1) · · ·Γ(αp1 + 1)

...

Bα,N
`,20 =

∑̀
p19 =0

p19∑
p18 =0

· · ·

p2∑
p1 =0

b`−p19N b
p19−p18
N · · ·b

p1
N

Γ(α` −αp19 + 1)Γ(αp19 −αp18 + 1) · · ·Γ(αp1 + 1)
.

(43)

One sees that the fractional differential equation under con-
sideration now takes the interesting form

c Dαy(x) +
ω

xα−β
c Dβy(x) + a0 + a1y(x) + a2y

2(x)

+ · · · + aN y
N (x) = 0, (44)

which clearly is a generalization of the fractional Lane-Emden
equation

c Dαy(x) +
ω

xα−β
c Dβy(x) + ym(x) = 0, (45)

in the sense that the new function f (y(x)) is a linear combination
of the classical polytropic term ym(x), 0 ≤ m ≤ N. Equation
(45) is the one considered in [28, Subsection 5.2.2] using the
power series method.

To make our calculations explicit and computationally in-
teresting we consider those elementary functions whose expan-
sion coefficients are explicitly known. Such elementary func-
tions to be examined with their explicit associated expansion
coefficients am, m = 0, 1, 2, 3, . . . , are the trigonometric, hyper-
bolic and exponential functions. These functions are enumer-
ated as follows.

(a) f (y(x)) = sin y(x), a2m = 0, a2m+1 = (−1)m/(2m + 1)!;

(b) f (y(x)) = cos y(x), a2m+1 = 0, a2m = (−1)m/(2m)!;

(c) f (y(x)) = sinh y(x), a2m = 0, a2m+1 = 1/(2m + 1)!;
269



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 270

(d) f (y(x)) = cosh y(x), a2m+1 = 0, a2m = 1/(2m)!;

(e) f (y(x)) = e±y(x), am = (±1)m/m!.

For extensive discussions of models involving these functions
as well as the physical interpretations of the solutions, see [12],
[13], [42].

We proceed with the computation of the approximate solu-
tion of the problem (10) according to whether the functions f (y)
are trigonometric, hyperbolic or exponential functions.

4.1. Lane-Emden Equation Involving Trigonometric Functions

In this subsection, we treat the equation (10) as well as
(34) for which the functions f (y(x)) are trigonometric func-
tions, namely, the sine function and cosine function. The Maclau-
rin series representations for these trigonometric functions are
employed.

4.1.1. The Special Case f (y(x)) = sin y(x)
In this case the fractional initial value problem (10) reduces

to a special one

c Dαy(x) +
ω

xα−β
c Dβy(x) + sin y(x) = 0,

y(0) = 1, y′(0) = 0.
(46)

In solving the initial value problem (46), we consider the series
expansion of sin y(x), namely,

sin y(x) =
N∑

m=0

a2m+1y
2m+1(x), a2m+1 =

(−1)m

(2m + 1)!
,

0 ≤ m ≤ N, N ∈ N0. (47)

In this case the problem (34) becomes

c Dαy(x) +
ω

xα−β
c Dβy(x) +

N∑
m=0

(−1)m
y2m+1(x)
(2m + 1)!

= 0

y(0) = 1, y′(0) = 0.

(48)

By Theorem 4.1, the analytical solution of the reduced frac-
tional problem (48) is then given by the Mittag-Leffler function
expansion

y(x) =
∞∑
`=0

B`N (α,β; ω)
Γ(α` + 1)

xα`, (49)

where the associated expansion coefficients are given explicitly
by

B
0
N (α,β; ω) := b

0
N = 1, (50)

B
`+1
N (α,β; ω) := b

`+1
N

= −
Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα

`,N

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
, ` ≥ 0. (51)

The coefficients Cα
`,N, ` ≥ 0, appearing in (51) take the formu-

lation

Cα`,N =
N∑

m=0

(−1)m

(2m + 1)!
Bα,N
`,2m+1. (52)

In order to see explicitly the values of the coefficients Cα
`,N ,

we consider the first values N ∈ N0 and this is presented as
follows.

(a) N = 0. In this case, equation (48) reduces to the initial
value problem

c Dαy(x) +
ω

xα−β
c Dβy(x) + y(x) = 0,

y(0) = 1, y′(0) = 0.
(53)

It follows from Theorem 4.1 that the solution of the initial
value problem (53) is given by the Mittag-Leffler expansion

y(x) =
∞∑
`=0

B`0(α,β; ω)
Γ(α` + 1)

xα`, (54)

where the expansion coefficients B`0 admit the explicit for-
mulation (` ≥ 0)

B
0
0(α,β; ω) = 1,

B
`+1
0 (α,β; ω) = −

Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα
`,0

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
, (55)

with the numbers Cα
`,0, ` ≥ 1, given by

Cα`,0 = B
α,0
`,1 = A

α
`,0 =

b`0
Γ(α` + 1)

. (56)

For further explicit calculations we proceed by assigning
special values to the parameters α,β and ω and this is done
in the following way as examples.

Example 4.1. Setting α = 2,β = 1,ω = 2, we see in this
case that equation (53) reduces to the initial value problem

y′′ +
2
x

y′ + y = 0, y(0) = 1, y′(0) = 0. (57)

It is seen from the series (54) that the solution of the initial
value problem (57) takes the formulation

y(x) =
∞∑
`=0

B`0(2, 1; 2)
Γ(2` + 1)

x2`, (58)

where the expansion coefficients B`0 = B
`
0(2, 1; 2) admit the

explicit formulation (` = 0, 1, 2, . . . )

B
0
0(2, 1; 2) = 1,

B
`+1
0 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,0

Γ(2` + 2) + 2Γ(2` + 1)
, (59)

270



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 271

and

C2`,0 = B
2,0
`,1 = A

2
`,0 =

b`0
Γ(2` + 1)

,

` = 0, 1, 2, . . . . (60)

Clearly, the first Mittag-Leffler expansion coefficients
B`0(2, 1; 2), ` = 1, 2, 3, 4 yield the following values.

B
1
0(2, 1; 2) = − C

2
0,0 ·

1
3

= −b00 ·
1
3

= −
1
3

B
2
0(2, 1; 2) = − C

2
1,0 ·

6
5

= −
B10

2!
·

6
5

=
1
6
·

6
5

=
1
5

B
3
0(2, 1; 2) = − C

2
2,0 ·

5!
7

= −
B20

4!
·

5!
7

= −
1
5!
·

5!
7

= −
1
7

B
4
0(2, 1; 2) = − C

2
3,0 · 560 = −

B30

6!
· 560 =

1
7!
· 560 =

1
9
.

(61)

Substituting these coefficients into (58) gives the solution of
the initial value problem (57):

y(x) = 1 −
1
3!

x2 +
1
5!

x4 −
1
7!

x6 +
1
9!

x8 + · · · =
sin x

x
.

(62)

Example 4.2. Indeed if we set α = 3/2,β = 1/2,ω = 2,
then one sees that equation (53) becomes the initial value
problem

c D
3
2 y(x) +

2
x

c D
1
2 y(x) + y(x) = 0,

y(0) = 1, y′(0) = 0.
(63)

One clearly sees from the series (54) that the solution of the
initial value problem (63) gives

y(x) =
∞∑
`=0

B`0

(
3
2 ,

1
2 ; 2

)
Γ
(

3`
2 + 1

) x 3`2 , (64)
where the expansion coefficients B`0 = B

`
0

(
3
2 ,

1
2 ; 2

)
admit

the Gamma function representation

B
0
0

(
3
2
,

1
2

; 2
)

= 1,

B
`+1
0

(
3
2
,

1
2

; 2
)

= −
Γ( 32 ` + 1)Γ(

3
2 ` + 2)C

3
2
`,0

Γ( 32 ` + 2) + 2Γ(
3
2 ` + 1)

, ` ≥ 0, (65)

and

C
3
2
`,0 = B

3
2 ,0
`,1 = B

3
2 ,0
`,1 =

b`0
Γ
(

3
2 ` + 1

), ` = 0, 1, 2, . . . . (66)

The first Mittag-Leffler expansion coefficients
B`0

(
3
2 ,

1
2 ; 2

)
, ` = 1, 2, 3, 4, 5, 6, are computed as follows:

B
1
0

(
3
2
,

1
2

; 2
)

= − C
3
2
0,0 ·

1
3

= −b00 ·
1
3

= −
1
3

B
2
0

(
3
2
,

1
2

; 2
)

= − C
3
2
1,0 ·

Γ( 52 ) · Γ(
7
2 )

Γ( 72 ) + 2Γ(
5
2 )

= −
B10

Γ( 52 )
5
√
π

12
=

5
27

B
3
0

(
3
2
,

1
2

; 2
)

= − C
3
2
2,0 ·

Γ(4) · Γ(5)
Γ(5) + 2Γ(4)

= −
B20

Γ(4)
· (4) = −

10
81

B
4
0

(
3
2
,

1
2

; 2
)

= − C
3
2
3,0 ·

Γ( 112 ) · Γ(
13
2 )

Γ( 132 ) + 2Γ(
11
2 )

= −
B30

Γ( 112 )
·

693
√
π

32

=
22
243

B
5
0

(
3
2
,

1
2

; 2
)

= − C
3
2
4,0 ·

Γ(7) · Γ(8)
Γ(8) + 2Γ(7)

= −
B40

Γ(7)
· (560)

= −
154
2187

B
6
0

(
3
2
,

1
2

; 2
)

= − C
3
2
5,0 ·

Γ
(

17
2

)
· Γ( 192 )

Γ( 192 ) + 2Γ(
17
2 )

= −
B50

Γ
(

17
2

) · 1640925√π
256

=
374
6561

. (67)

Upon substituting these coefficients into (64) yields the so-
lution

y(x) = 1 −
4x

3
2

9
√
π

+
5x3

162
−

64x9/2
15309

√
π

+
11x6

87480

−
512x15/2

57572775
√
π

+
187x9

1190427840
+ · · · . (68)

(b) N = 1. Indeed it is understood here that equation (48)
reduces to the initial value problem

c Dαy(x) +
ω

xα−β
c Dβy(x) + y(x) −

y3(x)
3!

= 0,

y(0) = 1, y′(0) = 0.
(69)

It is also seen here that the solution of the fractional prob-
lem (69) gives

y(x) =
∞∑
`=0

B`1(α,β; ω)
Γ(α` + 1)

xα`, (70)

with the expansion coefficients B`1 = B
`
1(α,β; ω) (` ≥ 0)

given by

B
0
1(α,β; ω) = 1,

B
`+1
1 (α,β; ω) = −

Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα
`,1

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
. (71)

Here the coefficients Cα
`,1, ` ≥ 0, on the right of (71) are

given by

Cα`,1 = B
α,1
`,1 −

1
3!

Bα,1
`,3 , (72)

271



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 272

where the numbers Bα,N
`,m , m = 1, 3, are as illustrated in (16)

and (27) respectively.

Example 4.3. In this case, another special Lane-Emden
equation can be obtained by setting α = 2,β = 1, ω = 2;
and as a result, equation (69) reduces to

y′′ +
2
x

y′ + y −
y3

3!
= 0, y(0) = 1, y′(0) = 0, (73)

whose solution admits the series representation

y(x) =
∞∑
`=0

B`1(2, 1; 2)
Γ(2` + 1)

x2`. (74)

From equations (71) and (72), we respectively have

B
0
1(2, 1; 2) = 1,

B
`+1
1 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,1

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0; (75)

and

C2`,1 =
b`1

Γ(2` + 1)

−
1
6

∑̀
p2 =0

p2∑
p1 =0

bp11 b
p2−p1
1 b

`−p2
1

Γ(2p1 + 1)Γ(2p2 − 2 p1 + 1)Γ(2` − 2p2 + 1)

 .
(76)

It is straightforward to see that the first Mittag-Leffler ex-
pansion coefficients B`1(2, 1; 2), ` = 1, 2, 3, 4, 5, 6, are cal-
culated as follows.

B
1
1(2, 1; 2) = −C

2
0,1 ·

1
3

= −
5b01
6
·

1
3

= −
5
6
·

1
3

= −
5

18

B
2
1(2, 1; 2) = −C

2
1,1 ·

6
5

= −
B11

3
·

6
5

=
5

54
·

6
5

=
1
9

B
3
1(2, 1; 2) = −C

2
2,1 ·

5!
7

= −

B2136 − B
1
1B

1
1

4

 · 1207
=

1
7776

·
120
7

=
5

2268

B
4
1(2, 1; 2) = −C

2
3,1 · 560 = −

 B311080 − B
1
1B

2
1

144

 · 560
= −

53
244944

· 560 = −
265
2187

B
5
1(2, 1; 2) = −C

2
4,1 ·

362880
11

= −

 B4160480 − B
1
1B

3
1

4320
−
B21B

2
1

3456

 · 36288011 = 287516038
B

6
1(2, 1; 2) = −C

2
5,1 ·

39916800
13

= −

 B515443200 − B
1
1B

4
1

241920
−
B21B

3
1

51840

 · 3991680013 = 2905085293.
(77)

Upon inserting the first coefficients (77) in (74) gives the
solution

y(x) =1 −
5

18
·

1
2!

x2 +
1
9
·

1
4!

x4 +
5

2268
·

1
6!

x6

−
265

2187
·

1
8!

x8 +
2875

16038
·

1
10!

x10 +
29050
85293

·
1

12!
x12 · · ·

=1 − 0.138889x2 + 0.00347222x4 + 0.0000030619x6

− 0.0000030052x8 + 0.0000000494x10

+ 0.0000000007x12 + · · · (78)

(c) N = 2. Another approximate solution of the problem(46)
can be obtained in this case in view of the reduced problem
(48):

c Dαy(x) +
ω

xα−β
c Dβy(x) + y(x) −

y3(x)
3!

+
y5(x)

5!
= 0,

y(0) = 1, y′(0) = 0.
(79)

Example 4.4. Similarly, in this case of N = 2, we spe-
cialise to the values α = 2,β = 1, ω = 2; to see that
equation (79) becomes

y′′ +
2
x

y′ + y−
y3

3!
+

y5

5!
= 0, y(0) = 1, y′(0) = 0.(80)

Clearly, the solution of the problem (80) takes the series
formula

y(x) =
∞∑
`=0

B`2(2, 1; 2)
Γ(2` + 1)

x2`, (81)

with the coefficients described as follows:

B
0
2(2, 1; 2) = 1,

B
`+1
2 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,2

Γ(2` + 2) + 2Γ(2` + 1)
,

C2`,2 = a1B
2,2
`,1 + a3B

2,2
`,3 + a5B

2,2
`,5 = B

2,2
`,1 −

1
3!

B2,2
`,3 +

1
5!

B2,2
`,5.

(82)

Here, one clearly understands that

B2,2
`,1 =

b`2
Γ(2` + 1)

,

B2,2
`,3 =

∑̀
p2 =0

p2∑
p1 =0

bp12 b
p2−p1
2 b

`−p2
2

Γ(2p1 + 1)Γ(2p2 − 2p1 + 1)Γ(2` − 2 p2 + 1)

B2,2
`,5 =

∑̀
p4 =0

p4∑
p3 =0

p3∑
p2 =0

p2∑
p1 =0

bp12 b
p2−p1
2 b

p3−p2
2 b

p4−p3
2 b

`−p4
2 b

`
2[Γ(2` − 2p4 + 1)]

−1

Γ(2p1 + 1)Γ(2p2 − 2p1 + 1)Γ(2 p3 − 2p2 + 1)Γ(2p4 − 2 p3 + 1)
.

(83)

272



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 273

As a result, the first Mittag-Leffler expansion coefficients
B`3(2, 1; 2), ` = 1, 2, 3, 4, are computed as follows.

B
1
2(2, 1; 2) = − C

2
0,2 ·

1
3

= −
101
120
·

1
3

= −
101
360

B
2
2(2, 1; 2) = − C

2
1,2 ·

6
5

= −
41B12
120

·
6
5

=
4141

36000

B
3
2(2, 1; 2) = − C

2
2,2 ·

5!
7

= −

 41B221440 − 19B
1
2B

1
2

480

 · 5!7
=

49591
18144000

B
4
2(2, 1; 2) = − C

2
3,2 · (560)

= −

 41B3243200 − 19B
1
2B

2
2

2880

 · 560 = − 16891139139968000
B

5
2(2, 1; 2) = − C

2
4,2 ·

362880
11

= −

 41B422419200 − 19B
1
2B

3
2

86400
−

19B22B
2
2

69120

 · 36288011
=

93349921751
513216000000

. (84)

Upon substituting these coefficients into the series formula
(81) gives the solution

y(x) =1 −
101
360
·

1
2!

x2 +
4141
36000

·
1
4!

x4 +
49591

18144000
·

1
6!

x6

−
16891139
139968000

·
1
8!

x8 +
93349921751

513216000000
·

1
10!

x10 + · · ·

=1 − 0.14027778x2 + 0.00479282x4

+ 0.000003796x6 + 0.000002993x8

+ 0.0000000501x10 + · · · (85)

4.1.2. The Special Case f (y(x)) = cos y(x)
It is seen here that the fractional Lane-Emden equation un-

der consideration is

c Dαy(x) +
ω

xα−β
c Dβy(x) + cos y(x) = 0,

y(0) = 1, y′(0) = 0.
(86)

Finding the approximate solution of the initial value prob-
lem (86) amounts to solving the problem

c Dαy(x) +
ω

xα−β
c Dβy(x) +

N∑
m=0

(−1)m

(2m)!
y2m(x) = 0,

y(0) = 1, y′(0) = 0,

(87)

where we have used the truncated power series

cos y(x) =
N∑

m=0

a2my
2m(x) =

N∑
m=0

(−1)m

(2m)!
y2m(x), N ∈ N0. (88)

Theorem 4.1 tells us that the analytical solution of the initial
value problem (87) admits the series solution

y(x) =
∞∑
`=0

B`N (α,β; ω)
Γ(α` + 1)

xα`, (89)

where we have

B
0
N (α,β; ω) := b

0
N = 1, (90)

B
`+1
N (α,β; ω) := b

`+1
N

= −
Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα

`,N

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
, ` ≥ 0. (91)

with

Cα`,N =
N∑

m=0

(−1)m

(2m)!
Bα,N
`,2m. (92)

We now proceed to the explicit computation of the first ex-
pansion coefficients B`N (α,β; ω). This is carried out by first con-
sidering the first values of N ∈ N0.

(a) N = 0. In this case, equation (86) reduces to the initial
value problem

c Dαy(x) +
ω

xα−β
c Dβy(x) + 1 = 0,

y(0) = 1, y′(0) = 0.
(93)

One sees from (89) that the solution of the initial value
problem (93) is given by

y(x) =
∞∑
`=0

B`0(α,β; ω)
Γ(α` + 1)

xα`, (94)

where (` ≥ 0)

B
0
0(α,β; ω) = 1,

B
`+1
0 (α,β; ω) = −

Γ(α` + 1)Γ(α(` + 1) −β + 1)Cα
`,0

Γ(α(` + 1) −β + 1) + ωΓ(α` + 1)
, (95)

with the numbers Cα
`,0, ` ≥ 0, given by

Cα`,0 = a0B
α,0
`,0 = B

α,0
`,0 = δ`0. (96)

Example 4.5. With the special values α = 2,β = 1 and ω =
2, we obtain a special case of the standard Lane-Emden
equation

y′′ +
2
x

y′ + 1 = 0, y(0) = 1, y′(0) = 0, (97)

whose solution is of the form

y(x) =
∞∑
`=0

B`0(2, 1; 2)
Γ(2` + 1)

x2`. (98)

Here we have

B
0
0(2, 1; 2) = 1,

B
`+1
0 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,0

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0, (99)

273



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 274

with C2
`,0 = δ`0. Interestingly one has the formula

B
1
0(2, 1; 2) =


−

1
3 , ` = 0,

0, ` = 1, 2, · · · ,
(100)

and as a result we obtain the closed form solution

y(x) = 1 −
1
6

x2. (101)

(b) N = 1. In this case, equation (87) reduces to the initial
value problem

c Dαy(x) +
ω

xα−β
c Dβy(x) + 1 −

y2

2!
= 0,

y(0) = 1, y′(0) = 0.
(102)

Example 4.6. Considering the the special values α = 2,β =
1 and ω = 2; equation (102) becomes the initial value prob-
lem

y′′ +
2
x

y′ + 1 −
y2

2!
= 0, y(0) = 1, y′(0) = 0, (103)

with the solution

y(x) =
∞∑
`=0

B`1(2, 1; 2)
Γ(2` + 1)

x2`. (104)

It follows that

B
0
1(2, 1; 2) = 1,

B
`+1
1 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,1

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0 (105)

with

C2`,1 =a0B
2,1
`,0 + a2B

2,1
`,2

=δ`0 −
1
2

∑̀
p1 =0

bp11 b
`−p1
1

Γ(2p1 + 1)Γ(2`− 2 p1 + 1)
.

Consequently, the first Mittag-Leffler expansion coefficients
B`1(2, 1; 2), ` = 1, 2, 3, 4, 5, are now computed as follows.

B
1
1(2, 1; 2) = − C

2
0,1 ·

1
3

= −
B01

2
·

1
3

= −
1
6

B
2
1(2, 1; 2) = − C

2
1,1 ·

6
5

=
B11

2
·

6
5

= −
1

10

B
3
1(2, 1; 2) = − C

2
2,1 ·

5!
7

=

B2124 + B
1
1B

1
1

8

 · 5!7 = − 184
B

4
1(2, 1; 2) = − C

2
3,1 · 560 =

 B31720 + B
1
1B

2
1

48

 · 560 = 527
B

5
1(2, 1; 2) = − C

2
4,1 ·

362880
11

=

 B4140320 + B
1
1B

3
1

1440
+
B21B

2
1

1152

 · 36288011 = 2960.
(106)

Upon using the coefficients (106), one obtains the solution

y(x) =1 −
1
6
·

1
2!

x2 −
1

10
·

1
4!

x4 −
1

84
·

1
6!

x6 +
5

27
·

1
8!

x8

+
29
60
·

1
10!

x10 + · · ·

=1 − 0.0833333x2 − 0.00416667x4 − 0.00001653x6

+ 0.00000459x8 + 0.00000013x10 + · · · (107)

(c) N = 2. Here, equation (87) becomes the initial value prob-
lem

c Dαy(x) +
ω

xα−β
c
Dβy(x) + 1 −

y2

2!
+

y4

4!
= 0,

y(0) = 1, y′(0) = 0.
(108)

Example 4.7. Consider the problem

y′′+
2
x

y′+1−
y2

2!
+

y4

4!
= 0, y(0) = 1, y′(0) = 0.(109)

The solution of (109) takes the form

y(x) =
∞∑
`=0

B`2(2, 1; 2)
Γ(2` + 1)

x2`, (110)

where the coefficients B`2(2, 1; 2) are given by

B
0
2(2, 1; 2) = 1,

B
`+1
2 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,2

Γ(2` + 2) + 2Γ(2` + 1)
, (` ≥ 0)

C2`,2 = a0B
2,2
`,0 + a2B

2,2
`,2 + a4B

2,2
`,4 = B

2,2
`,0 −

1
2!

B2,2
`,2 +

1
4!

B2,2
`,4,

(111)

with

B2,2
`,0 = δ`0, B

2,2
`,2 =

∑̀
p1 =0

bp12 b
`−p1
2

Γ(2p1 + 1)Γ(2`− 2 p1 + 1)

B2,2
`,4 =

∑̀
p3 =0

p3∑
p2 =0

p2∑
p1 =0

bp12 b
p2−p1
2 b

p3−p2
2 b

`−p3
2

Γ(2p1 + 1)Γ(2p2 − 2p1 + 1)Γ(2 p3 − 2p2 + 1)Γ(2`− 2p3 + 1)
.

(112)

The first Mittag-Leffler expansion coefficients B`+12 (2, 1; 2), ` =
0, 1, 2, . . . are computed as follows.

B
1
2(2, 1; 2) = −C

2
0,2 ·

1
3

= −

1 − B022 + B
0
2

24

 · 13 = −1372
B

2
2(2, 1; 2) = −C

2
1,2 ·

6
5

= −

−B122 + B
1
2

24

 · 65 = − 1431440
274



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 275

B
3
2(2, 1; 2) = − C

2
2,2 ·

5!
7

=

 11B22288 + 11B
1
2B

1
2

96

 · 1207
= −

143
145152

B
4
2(2, 1; 2) = − C

2
3,2 · 560 =

 11B328640 + 11B
1
2B

2
2

576

 · 560
=

26741
139968

B
5
2(2, 1; 2) = − C

2
4,2 ·

362880
11

=

 11B42483840 + 11B
1
2B

3
2

17280
+

11B22B
2
2

1152

 · 36288011
=

6059911
14929920

. (113)

Upon substituting these coefficients in (110) yields the se-
ries solution

y(x) =1 −
13
72
·

1
2!

x2 −
143
1440

·
1
4!

x4 −
143

145152
·

1
6!

x6

+
26741

139968
·

1
8!

x8 +
6059911
14929920

·
1

10!
x10 + · · ·

=1 − 0.09027778x2 − 0.00413773x4 − 0.00000137x6

+ 0.000004738x8 + 0.000000112x10 + · · · . (114)

4.2. Lane-Emden Equation Involving Hyperbolic Functions
This subsection discusses the approximate solution of the

initial value problems (10) and (34) for which the functions
f (y(x)) are the hyperbolic sine function and the hyperbolic co-
sine function. We also make use of the Maclaurin series repre-
sentations for these hyperbolic functions.

4.2.1. The Special Case f (y(x)) = sinh y(x)
We consider the fractional Lane-Emden initial value prob-

lem

c Dαy(x) +
ω

xα−β
c
Dβy(x) + sinh y(x) = 0,

y(0) = 1, y′(0) = 0.
(115)

In solving the initial value problem (115), we are interested in
the series expansion of sinh y(x), namely (0 ≤ m ≤ N, N ∈ N0),

sinh y(x) =
N∑

m=0

a2m+1y
2m+1(x), a2m+1 =

1
(2m + 1)!

, (116)

so that finding the approximate solution of the problem (115)
amounts to solving the corresponding problem

c Dαy(x) +
ω

xα−β
c
Dβy(x) +

N∑
m=0

1
(2m + 1)!

y2m+1(x) = 0,

y(0) = 1, y′(0) = 0.

(117)

We proceed to obtaining the solutions of the initial value
problem (117) for the first values of N ∈ N.

(a) N = 0. Clearly, in this case a1 = 1 and as a result the
resulting initial value problem coincides with the problem
in (53).

(b) N = 1. Here, we have from (117) the problem

c Dαy(x) +
ω

xα−β
c
Dβy(x) + y(x) +

y3(x)
3!

= 0,

y(0) = 1, y′(0) = 0.
(118)

Example 4.8. Consider the initial value problem

y′′ +
2
x

y′ + y +
y3

3!
= 0, y(0) = 1, y′(0) = 0 (119)

Clearly the expansion coefficients appearing in the series
solution of (119) are given by

B
0
1(2, 1; 2) = 1,

B
`+1
1 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,1

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0. (120)

with

C2`,1 = a1B
2,1
`,1 + a3B

2,1
`,3

=
b`1

Γ(2` + 1)

+
1
6

∑̀
p2 =0

p2∑
p1 =0

bp11 b
p2−p1
1 b

`−p2
1

Γ(2p1 + 1)Γ(2p2 − 2 p1 + 1)Γ(2` − 2p2 + 1)

 .
(121)

The first Mittag-Leffler expansion coefficients B`1(2, 1; 2), ` =
0, 1, 2, 3, 4, 5, are presented as follows.

B
1
1(2, 1; 2) = −C

2
0,1 ·

1
3

= −
7B01

6
·

1
3

= −
7

18

B
2
1(2, 1; 2) = −C

2
1,1 ·

6
5

= −
2B11

3
·

6
5

=
14
45

B
3
1(2, 1; 2) = −C

2
2,1 ·

5!
7

= −

B2118 − B
1
1B

1
1

24

 · 1207 = −131324
B

4
1(2, 1; 2) = −C

2
3,1 · 560 = −

 B31540 + B
1
1B

2
1

144

 · 560 = 19462187
B

5
1(2, 1; 2) = −C

2
4,1 ·

362880
11

= −

 B4130240 + B
1
1B

3
1

4320
+
B21B

2
1

3456

 · 36288011 = −24821380190
(122)

Upon substituting the first coefficients (122) we obtain the

275



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 276

solution

y(x) =
∞∑
`=0

B`1(2, 1; 2)
Γ(2` + 1)

x2`

=1 −
7

18
·

1
2!

x2 +
14
45
·

1
4!

x4 −
131
324
·

1
6!

x6 +
1946
2187

·
1
8!

x8

−
248213
80190

·
1

10!
x10 + · · ·

=1 − 0.19444444x2 + 0.01296296x4 − 0.00056156x6

+ 0.000022069x8 − 0.000000853x10 + · · · . (123)

(c) N = 2. Indeed it is easy to understand here that from (117)
one obtains the initial value problem

c Dαy(x) +
ω

xα−β
c
Dβy(x) + y(x) +

y3(x)
3!

+
y5(x)

5!
= 0,

y(0) = 1, y′(0) = 0.
(124)

Example 4.9. Given the problem

y′′+
2
x

y′+y+
y3

3!
+

y5

5!
= 0, y(0) = 1, y′(0) = 0.(125)

Clearly, the expansion coefficients of the series solution of
(125) are

B
0
2(2, 1; 2) = 1,

B
`+1
2 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,2

Γ(2` + 2) + 2Γ(2` + 1)
,

C2`,2 = a1B
2,2
`,1 + a3B

2,2
`,3 + a5B

2,2
`,5 = B

2,2
`,1 +

1
3!

B2,2
`,3 +

1
5!

B2,2
`,5,

(126)

where

B2,2
`,1 =

b`2
Γ(2` + 1)

,

B2,2
`,3 =

∑̀
p2 =0

p2∑
p1 =0

bp12 b
p2−p1
2 b

`−p2
2

Γ(2p1 + 1)Γ(2p2 − 2 p1 + 1)Γ(2` − 2p2 + 1)

B2,2
`,5 =

∑̀
p4 =0

p4∑
p3 =0

p3∑
p2 =0

p2∑
p1 =0

bp12 b
p2−p1
2 b

p3−p2
2 b

p4−p3
2 b

`−p4
2 [Γ(2` − 2 p4 + 1)]

−1

Γ(2p1 + 1)Γ(2p2 − 2 p1 + 1)Γ(2p3 − 2p2 + 1)Γ(2p4 − 2p3 + 1)
.

(127)

The first expansion coefficients B`1(2, 1; 2), ` = 1, 2, 3, 4, 5

are computed as follows.

B
1
2(2, 1; 2) = −C

2
0,2 ·

1
3

= −
47b02
40
·

1
3

= −
47

120

B
2
2(2, 1; 2) = −C

2
1,2 ·

6
5

= −
27B12

3
·

6
5

=
1269
4000

B
3
2(2, 1; 2) = −C

2
2,2 ·

5!
7

= −

 9B22160 − 7B
1
2B

1
2

160

 · 1207
= −

282893
672000

B
4
2(2, 1; 2) = −C

2
3,2 · 560 = −

 3B321600 + 7B
1
2B

2
2

960

 · 560
=

4747
5000

B
5
2(2, 1; 2) = −C

2
4,2 ·

362880
11

= −

 3B4289600 + 7B
1
2B

3
2

28800
+

7B22B
2
2

23040

 · 36288011 = −2379140611704000000
(128)

Upon inserting the first coefficients into the series solution
one obtains

y(x) =
∞∑
`=0

B`2(2, 1; 2)
Γ(2` + 1)

x2`

= 1−
47

120
·

1
2!

x2 +
1269
4000

·
1
4!

x4−
282893
672000

·
1
6!

x6 +
4747
5000

·
1
8!

x8

−
2379140611
704000000

·
1

10!
x10 + · · ·

= 1 − 0.19583333x2 + 0.01321875x4 − 0.00058468x6

+ 0.000023547x8 − 0.000000931x10 + · · · . (129)

4.2.2. The Special Case f (y(x)) = cosh y(x)
We consider the fractional Lane-Emden equation

c Dαy(x) +
ω

xα−β
c
Dβy(x) + cosh y(x) = 0,

y(0) = 1, y′(0) = 0.
(130)

In order to solve the initial value problem (130), we consider
the series expansion of cosh y(x):

cosh y(x) =
N∑

m=0

a2my
2m(x) =

N∑
m=0

y2m(x)
(2m)!

, 0 ≤ m ≤ N, N ∈ N0,

(131)

and as a result we have the associated problem

c Dαy(x) +
ω

xα−β
c
Dβy(x) +

N∑
m=0

y2m(x)
(2m)!

= 0,

y(0) = 1, y′(0) = 0.

(132)

For N = 0, we observe that the resulting initial value prob-
lem coincides with the problem in (93).

276



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 277

Example 4.10. Given the initial value problem

y′′ +
2
x

y′ + 1 +
y2

2!
= 0, y(0) = 1, y′(0) = 0. (133)

It follows immediately that the solution to the problem (133) is
given by

y(x) =1 −
1
2
·

1
2!

x2 +
3
10
·

1
4!

x4 −
3
4
·

1
6!

x6 +
7
3
·

1
8!

x8

−
2877
220

·
1

10!
x10 + · · ·

=1 − 0.25000000x2 + 0.01250000x4 − 0.00104167x6

+ 0.00005787x8 − 0.00000360x10 + · · · (134)

Example 4.11. Consider the problem

y′′ +
2
x

y′ + 1 +
y2

2!
+

y4

4!
= 0, y(0) = 1, y′(0) = 0. (135)

One easily sees that the solution to the problem (135) gives

y(x) =1 −
37
72
·

1
2!

x2 +
481

1440
·

1
4!

x4 −
126503
145152

·
1
6!

x6

+
406445
139968

·
1
8!

x8 −
256054097
14929920

·
1

10!
x10 · · ·

=1 − 0.25694444x2 + 0.01391782x4 − 0.00121045x6

+ 0.00007202x8 − 0.00000473x10 + · · · . (136)

4.3. The Case of Exponential Functions f (y) = exp(±y)
In this subsection, we compute the approximate solution of

the fractional Lane-Emden problem (10) as well as the initial
problem (34) for which the functions f (y(x)) are exponential
functions. The Maclaurin series expansion for these exponen-
tial functions are applied.

4.3.1. Fractional Isothermal Gas Sphere Equation
The fractional isothermal gas sphere equation is given by

c Dαy(x) +
ω

xα−β
c Dβy(x) + ey(x) = 0,

y(0) = 0, y′(0) = 0.
(137)

Also in this case, we look for the approximate solution of the
problem (137) by solving the complementary equation

c Dαy(x) +
ω

xα−β
c Dβy(x) +

N∑
m=0

ym(x)
m!

= 0,

y(0) = 0, y′(0) = 0.

(138)

The problem

y′′ +
2
x

y′ + ey = 0, y(0) = 0, y′(0) = 0 (139)

models the isothermal gas spheres where the temperature re-
mains constant and the index m is infinite ([42]).

Example 4.12. Given the problem

y′′ +
2
x

y′ + 1 = 0, y(0) = 0, y′(0) = 0. (140)

It is straightforward to see by the initial conditions in (140)
that the expansion coefficients appearing in the series solution
of (140) are given by

B
0
0(2, 1; 2) = 0,

B
`+1
0 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,0

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0, (141)

with C2
`,0 = δ`0. Thus one has the Mittag-Leffler expansion co-

efficient

B
1
0(2, 1; 2) =


−

1
3 , ` = 0,

0, ` = 1, 2, · · · .
(142)

Consequently we obtain the closed form solution

y(x) =
∞∑
`=0

B`0(2, 1; 2)
Γ(2` + 1)

x2` = −
1
6

x2. (143)

Example 4.13. N = 1: It is seen here that

y′′ +
2
x

y′ + 1 + y = 0, y(0) = 0, y′(0) = 0. (144)

The coefficients of the series solution of (144) has the formula-
tion

B
0
1(2, 1; 2) = 0,

B
`+1
1 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,1

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0 (145)

with

C2`,1 = δ`0 +
b`1

Γ(2` + 1)
. (146)

Further simplification gives the following first expansion coef-
ficients B`1(2, 1; 2), ` = 1, 2, 3, 4, 5.

B
1
1(2, 1; 2) = − C

2
0,1 ·

1
3

= −(1 + B01) ·
1
3

= −
1
3

B
2
1(2, 1; 2) = − C

2
1,1 ·

6
5

= −
B11

2
·

6
5

=
1
5

B
3
1(2, 1; 2) = − C

2
2,1 ·

5!
7

= −
B21

24
·

120
7

= −
1
7

B
4
1(2, 1; 2) = − C

2
3,1 · 560 = −

B31

720
· 560 =

1
9

B
5
1(2, 1; 2) = − C

2
4,1 ·

362880
11

= −
B41

40320
·

362880
11

= −
1
11
(147)

277



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 278

Upon substituting we obtain the solution

y(x) =
∞∑
`=0

B`1(2, 1; 2)
Γ(2` + 1)

x2`

= −
1
3
·

1
2!

x2 +
1
5
·

1
4!

x4 −
1
7
·

1
6!

x6

+
1
9
·

1
8!

x8 −
1
11
·

1
10!

x10 + · · ·

= −
1
3!

x2 +
1
5!

x4 −
1
7!

x6 +
1
9!

x8 −
1

11!
x10 + · · · (148)

Example 4.14. N = 2 : Consider the initial value problem

y′′ +
2
x

y′ + 1 + y +
y2

2!
= 0, y(0) = 0, y′(0) = 0. (149)

The series solution of (149) has the following expansion coeffi-
cients.

B
0
2(2, 1; 2) = 0,

B
`+1
2 (2, 1; 2) = −

Γ(2` + 1)Γ(2` + 2)C2
`,2

Γ(2` + 2) + 2Γ(2` + 1)
, ` ≥ 0, (150)

where

C2`,2 = δ`0 +
b`2

Γ(2` + 1)

+
1
2

∑̀
p1 =0

bp12 b
`−p1
2

Γ(2p1 + 1)Γ(2`− 2 p1 + 1)

 . (151)
Simplifying further gives the following first expansion coeffi-
cients B`2(2, 1; 2), ` = 0, 1, 2, . . . .

B
1
2(2, 1; 2) = −C

2
0,2 ·

1
3

= −1 ·
1
3

= −
1
3

B
2
2(2, 1; 2) = −C

2
1,2 ·

6
5

= −
B12

2
·

6
5

=
1
5

B
3
2(2, 1; 2) = −C

2
2,2 ·

5!
7

= −

B2224 + B
1
2B

1
2

8

 · 1207 = − 821
B

4
2(2, 1; 2) = −C

2
3,2 · 560 = −

 B32720 + B
1
2B

2
2

48

 · 560 = 2927
B

5
2(2, 1; 2) = −C

2
4,2 ·

362880
11

= −

 B4240320 + B
1
2B

3
2

1440
+
B22B

2
2

576

 · 36288011 = −403208173 . (152)
As a result we obtain the solution

y(x) =
∞∑
`=0

B`2(2, 1; 2)
Γ(2` + 1)

x2`

= −
1
3
·

1
2!

x2 +
1
5
·

1
4!

x4−
8
21
·

1
6!

x6 +
29
27
·

1
8!

x8−
40320
8173

·
1

10!
x10 +· · ·

= −
1
3!

x2 +
1
5!

x4 −
8

3 · 7!
x6 +

29
3 · 9!

x8 −
40320

743 · 11!
x10 + · · ·

(153)

4.3.2. Richardson’s Model of Thermionic Current
Here, we consider the initial value problem

c Dαy(x) +
ω

xα−β
c
Dβy(x) + e−y(x) = 0,

y(0) = 0, y′(0) = 0.
(154)

The problem (Richardson’s model of thermionic current)

y′′ +
2
x

y′ + e−y = 0, y(0) = 0, y′(0) = 0 (155)

models the density and electric force of an electron gas in the
neighbourhood of a hot body in thermal equilibrium ( [13],
[42]).

Clearly, for N = 0 with α = 2,β = 1 and ω = 2, (154)
coincides with (140).

Example 4.15. Consider the problem

y′′ +
2
x

y′ + 1 − y = 0, y(0) = 0, y′(0) = 0. (156)

As in the preceding calculations we easily see that the initial
value problem (156) has the solution

y(x) = −
1
3
·

1
2!

x2−
1
5
·

1
4!

x4−
1
7
·

1
6!

x6−
1
9
·

1
8!

x8−
1
11
·

1
10!

x10+· · ·

= −
1
3!

x2 −
1
5!

x4 −
1
7!

x6 −
1
9!

x8 −
1

11!
x10 + · · · . (157)

Example 4.16. N = 2 : Given the initial value problem

y′′ +
2
x

y′ + 1 − y +
y2

2!
= 0, y(0) = 0, y′(0) = 0. (158)

It is also easy to show that the solution to (158) yields

y(x) = −
1
3
·

1
2!

x2 −
1
5
·

1
4!

x4 −
8
21
·

1
6!

x6 −
29
27
·

1
8!

x8

−
40320
8173

·
1

10!
x10 + · · ·

= −
1
3!

x2 −
1
5!

x4 −
8

3 · 7!
x6 −

29
3 · 9!

x8 −
40320

743 · 11!
x10 + · · · .

(159)

5. Concluding Remarks

In this paper, we have used the Mittag-Leffler function ex-
pansion method to find approximate solutions of a class of frac-
tional Lane-Emden equation whose nonlinear forms of f (y) are
expressible as f (y(x)) = a0 + a1y(x) + a2y2(x) + · · ·+ amym(x) +
· · · + aN yN (x), 0 ≤ m ≤ N, N ∈ N0; the values of the expansion
coefficients am, 0 ≤ m ≤ N, were explicitly provided. We con-
sidered the cases for which the functions f (y(x)) are trigono-
metric, hyperbolic and exponential functions. In all these cases,
the Lane-Emden equations for the special values N = 0; α =
2,β = 1,ω = 2, coincide with the classical and standard ones
((57), (97), (140)) and consequently the corresponding solu-
tions are given in closed forms ((62), (101), (143)).

In the case where the functions f (y) are exponential func-
tions, our results can be compared with those of [11, equations

278



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 279

(59) and (70)] (see also [25, equation (25)]). The other non-
linear forms of f (y) are the trigonometric and hyperbolic func-
tions. Our approximate solutions in the case of trigonometric
functions are comparable with the solutions [11, equations (81)
and (83) ]; and in the case of hyperbolic functions, see [11,
equations (89) and (91)] (see also [22]).

The method employed in this paper can be applied to other
similar cases in applied sciences where the models are given
as strongly nonlinear ordinary differential equations. By exten-
sion, the method, can therefore, be accurately and reliably used
to compute approximate solutions of nonlinear fractional differ-
ential equations of Lane-Emden type where the nonlinear forms
of f (y) involve several other special functions.

Acknowledgements

The authors thank the anonymous referees for their careful
reading of the manuscript and valuable comments for the im-
provement of the quality of the paper.

References

[1] R. Garrappa, E. Kaslik & M. Popolizio, “Evaluation of fractional
integrals and derivatives of elementary functions: overview and
tutorial”, Mathematics 7, (2019) 407.

[2] G. Jumarie, “A Fokker-Planck equation of fractional order with
respect to time”, Journal of Math. Physics 33 (1992) 3536.

[3] G. Jumarie, “Fractional Fokker-Planck equation, solutions and
applications”, Physical Review, 63 (2001) 1.

[4] G. Jumarie, “Schrödinger equation for quantum-fractal space-
time of order n via the complex-valued fractional Brownian mo-
tion”, Intern. Y. of Modern Physics A 16 (2001) 5061.

[5] A. A. Kilbas, Theory and applications of fractional differential
equations, Elsevier, (2006).

[6] F. Mainardi & R. Gorenflo, “Time-fractional derivatives in re-
laxation processes: A tutorial survey”, Fractional Calculus and
Applied Analysis 10 (2007) 269.

[7] M. D. Ortigueira, Fractional calculus for scientists and engineers
Springer (2011).

[8] I. Podlubny, Fractional differential equations, Academic Press,
San Diego, (1999).

[9] M. Žecová & J. Terpák, “Heat conduction modeling by using
fractional-order derivatives”, Applied Mathematics and Compu-
tation 257 (2015) 365.

[10] K. Diethelm, The analysis of fractional differential equations:
an application-oriented exposition using differential operators of
Caputo type, Springer-Verlag Berlin Heidelberg (2010).

[11] A. M. Wazwaz, “A new algorithm for solving differential equa-
tions of Lane-Emden type”, Appl. Math. Comput. 118 (2001)
287.

[12] S. Chandrasekhar, Introduction to the Study of Stellar Structure,
Dover, New York, (1967).

[13] O. U. Richardson, The Emission of Electricity from Hot Bodies,
Longman, Green and Co., London, New York, (1921).

[14] C. Mohan & A. R. Al-Bayaty, “Power series solutions of the
Lane-Emden equation”, Astrophysics and Space Science 73
(1980) 227.

[15] J. I. Ramos, “Series approach to the Lane-Emden equation and
comparison with the homotopy perturbation method”, Chaos,
Solitons & Fractals 38 (2008) 400.

[16] S. K. Vanani & A. Aminataei, “On the numerical solutions of dif-
ferential equations of Lane-Emden type”, Computers and Mathe-
matics with Applications 59 (2010) 2815.

[17] M. O. Ogunniran, O. A. Tayo, Y. Haruna & A. F. Adebisi, “Lin-
ear stability analysis of Runge-Kutta methods for singular Lane-
Emden equations”, Journal of the Nigerian Society of Physical
Sciences 2 (2020) 134.

[18] E. O. Adeyefa & O. S. Esan, “Exponentially fitted Chebyshev
based algorithm as second order initial value solver”, Journal of
the Nigerian Society of Physical Sciences 2 (2020) 51.

[19] S. E. Fadugba, S.N. Ogunyebi & B.O. Falodun, “An examina-
tion of a second order numerical method for solving initial value
problems”, Journal of the Nigerian Society of Physical Sciences
2 (2020) 120.

[20] M. S. Mechee & N. Senu, “Numerical study of fractional differ-
ential equations of Lane-Emden type by method of collocation”,
Applied Mathematics 3 (2012) 851.

[21] A. Akgül, M. Inc, E. Karatas & D. Baleanu, “Numerical solu-
tions of fractional differential equations of Lane-Emden type by
an accurate technique”, Advances in Difference Equations 2015
(2015) 220.

[22] A. Saadatmandi, A. Ghasemi-Nasrabady & A. Eftekhari, “Nu-
merical study of singular fractional Lane-Emden type equations
arising in astrophysics”, J. Astrophys. Astr. (2019) 40 27.

[23] P. K. Sahu & B. Mallick, “Approximate solution of fractional
order Lane–Emden type differential equation by orthonormal
Bernoulli’s polynomials”, Int. J. Appl. Comput. Math (2019) 5
89.

[24] M. I. Nouh & E. A.-B. Abdel-Salam, “Approximate Solution to
the Fractional Lane–Emden Type Equations”, Iran J. Sci. Tech.
Trans. Sci. 42, (2018) 2199.

[25] A. M. Malik & O.H. Mohammed, “Two efficient methods for
solving fractional Lane–Emden equations with conformable frac-
tional derivative”, Journal of the Egyptian Mathematical Society
(2020) 28.

[26] B. Cǎruntu, C. Bota, M. Lǎpǎdat & M. S. Paşca, “Polyno-
mial least squares method for fractional Lane–Emden equations”,
Symmetry 11 (2019) 479.

[27] J. Davila, L. Dupaigne & J. Wei, “On the fractional Lane-Emden
equation”, Trans. Am. Math. Soc. 369 (2017) 6087.

[28] C. Milici, G. Drǎgǎnescu & J.T. Machado, “Introduction to frac-
tional differential equations”, Nonlinear Systems and Complex-
ity, 25 (2019).

[29] U. Saeed, “Haar Adomian method for the solution of fractional
nonlinear Lane-Emden type equations arising in astrophysics”,
Taiwanese Journal of Mathematics 21 (2017) 1175.

[30] O. A. Uwaheren, A. F. Adebisi & O. A. Taiwo, “Perturbed col-
location method for solving singular multi-order fractional dif-
ferential equations of Lane-Emden type”, Journal of the Nigerian
Society of Physical Sciences 2 (2020) 141.

[31] A. A. M. Arafa & S. Z. Rida, A. A. Mohammadein, H. M. Ali,
“Solving nonlinear fractional differential equation by generalized
Mittag-Leffler function method”, Communications in Theoretical
Physics 59 (2013) 661.

[32] A. Atangana & A. Secer, “A note on fractional order derivatives
and table of fractional derivatives of some special functions”, Ab-
stract and Applied Analysis 2013 (2013) 1.

[33] M. Caputo & M. Fabrizio, “A new definition of fractional deriva-
tive without singular kernel”, Progr. Fract. Differ. Appl. 1 (2015)
73.

[34] S. Das, “Functional Fractional Calculus”, Springer (2011).
[35] M. Davison & C. Essex, “Fractional differential equations and

279



Awonusika et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 265–280 280

initial value problems”, The Mathematical Scientist 23 (1998)
108.

[36] S. E. Fadugba, “Solution of fractional order equations in the do-
main of the Mellin transform”, Journal of the Nigerian Society of
Physical Sciences 1 (2019) 138.

[37] R. Herrmann, Fractional calculus, an introduction for Phycists,
World Scientific, (2011).

[38] J. Jumarie, “Modified Riemann-Liouville derivative and frac-
tional Taylor series of nondifferentiable functions further results”,
Computers and Mathematics with Applications 51 (2006) 1367.

[39] J. Jumarie, “Table of some basic fractional calculus formulae

derived from a modified Riemann-Liouville derivative for non-
differentiable functions”, Applied Mathematics Letters 22 (2009)
378.

[40] E. C. de Oliveira & J. A. T. Machado, “A review of definitions
for fractional derivatives and integral”, Mathematical Problems
in Engineering 2014 (2014) 1.

[41] F. Mainardi, Fractional calculus and waves in linear viscoelas-
ticity, Imperial College Press (2010).

[42] H. T. Davis, Introduction to Nonlinear Differential and Integral
Equations, Dover, New York, (1962).

280