J. Nig. Soc. Phys. Sci. 5 (2023) 1075

Journal of the
Nigerian Society

of Physical
Sciences

Collocation Method for the Numerical Solution of Multi-Order
Fractional Differential Equations

G. Ajileyea, A. A. Jamesb,∗

aDepartment of Mathematics and Statistics, Federal University Wukari, Taraba State, Nigeria.
b Department of Mathematics and Statistics, American University of Nigeria, Yola, Adamawa State, Nigeria.

Abstract

This study presents a collocation approach for the numerical integration of multi-order fractional differential equations with initial conditions in
the Caputo sense. The problem was transformed from its integral form into a system of linear algebraic equations. Using matrix inversion, the
algebraic equations are solved and their solutions are substituted into the approximate equation to give the numerical results. The effectiveness
and precision of the method were illustrated with the use of numerical examples.

DOI:10.46481/jnsps.2023.1075

Keywords: Keywords: Differential equation, Fractional derivatives, Approximate solution, Power series.

Article History :
Received: 17 September 2022
Received in revised form: 24 May 2023
Accepted for publication: 25 May 2023
Published: 11 June 2023

© 2023 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license

(https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

Communicated by: Tolulope Latunde

1. Introduction

In the fields of mathematics, physics, chemistry, and
engineering, differential and integral equations involving
fractions are of the utmost significance. The use of functional
equations, such as ordinary and partial differential equations,
is typical when applying mathematics to the modeling of
problems arising in the real world. In the early 1900s, Italian
mathematician Vito Volterra came up with a whole new sort
of equation that came to be known as integro-differential
equations in order to investigate the phenomenon of population
expansion. In these types of equations, one or more derivatives
of the function whose value is unknown is placed under the
integral sign. Integro-differential equations can be found in a

∗Corresponding author tel. no: +2348077092831
Email address: adewalejames974@gmail.com (A. A. James )

variety of mathematical formulations of physical phenomena.
Additionally, these equations can be found in the modeling of
certain phenomena in the fields of science and engineering.
For instance, the equations of kinetics that support the kinetic
theory of rarefied gases, plasma, radiation transmission, and
coagulation are some examples. [1]. Some of the numerical
solution of fractional differential equations developed in the
literature include: Perturbed collocation method [2], Adomian
decompositions method by [3-5], Collocation method by [6-9],
Chebyshev- Gelerkin method [10], Bernoulli matrix method
[11], Differential transform method [12], Pseudospectral
method [13], Bernstein Polynomials method [14, 15], the
Mellin transform approach [16]. [17] utilized a numerical
approach based on the boubaker polynomial to generate
approximate numerical solutions to the multi-order fractional
differential equations. Their decision was to use an operational
matrix for fractional integration based on boubakar polynomi-

1



Ajileye & James / J. Nig. Soc. Phys. Sci. 5 (2023) 1075 2

als. Collocation approach for the computational solution of
fredholm-volterra fractional order of integro-differential equa-
tions was presented by [18]. They solved the problem by first
obtaining the linear integral form of it and then transforming
it into a system of linear algebraic equations by making use of
conventional collocation points.

In this research, the collocation method is utilized to solve
multi-order fractional differential equations of the form

Dβy(x) =
N∑

j=0

q j(x)D
α j y(x) + h(x) (1)

subject to the initial condition

y( j)(a j) = λ j, j = 0, 1, ..., n − 1, n ∈ N β > αN, (2)

where y(x) is the unknown function, Dα j and Dβ are the Ca-
puto’s derivative, h(x) is the force known -prior. q j(x) is the
known function, a j and λ j are known constants.

2. Basic Definitions

In this section, we present certain definitions and fundamental
ideas of fractional calculus for the purpose of the formulation
of the problem that has been presented.
Definition 2.1: The Caputo derivative with order α > 0 of the
given function f (x), x ∈ (a, b) is defined as

C
x D

α
a y(x) =

1
Γ(m −α)

∫ x
a

(x − s)m−α−1y(m)(s)d s (3)

where m − 1 ≤ α ≤ m, m ∈ N, x > 0
Definition 2.2: Let (an) , n ≥ 0 be a sequence of real numbers.
The power series in x with coefficients an is an expression

y(x) = a0+a1 x+a2 x
2
+a3 x

3
+· · ·aN x

N
=

N∑
n=0

an x
n

= φ(x) A(4)

where φ(x) = [1 x x2 · · · xN ], A = [a0 a1 · · · aN ]T
then y(x, n) = xnA, n = 0(1)N, n ∈ Z+

Definition 2.3: Standard Collocation Method (SCM). This
method is used to determine the desired collocation points
within an interval. i.e [a,b] and is given by

xi = a +
(b − a)i

N
, i = 1, 2, 3, ...N (5)

Definition 2.4: Let y(x) be a continuous function, then

0 I
β
x

(
C
0 D

β
xy(x)

)
= y(x) −

N∑
k=0

y(k)(0)
k!

xk (6)

where m − 1 < β < 1
Definition 2.5: Let p(s) be an integrable function, then

0 I
β
x ( p(s)) =

1
Γ(β)

∫ x
0

(x − s)β−1 p(s)d s (7)

Definition 2.6: The Riemann -Liouville derivative of order α >
0 with n − 1 < α < n of the power function f (t) = t p−α is given
by

Dαt p =
Γ( p + 1)

Γ( p −α + 1)
t p−α (8)

3. Mathematical Background

In this part, we create a collocation approach for numeri-
cally solving multi-order fractional differential equations utiliz-
ing power series polynomials as the basis function.
Lemma (3.1) (Integral form)
Let y(x) be a solution to (1) subject to (2), the integral form is

y(x) = W(x) +
N∑

j=0

1
Γ(m j −α j)

1
Γ(β)

(9)

×

∫ x
0

(x − s)β−1q j(s)
[∫ s

0
(s − t)m j−α j−1y(m j )(t)dt

]
d s

where

W(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(β)

∫ x
0

(x − s)β−1h(s)d s

Proof. Multiply equation (1) by 0 I
β
x (.) gives

0 I
β
x

(
Dβy(x)

)
= 0 I

β
x

 N∑
j=0

q j(x)D
α j y(x)

+ 0 Iβx (h(x)) (10)
using (6) on equation (10) gives

y(x) =
N∑

k=0

y(k)(0)
k!

xk + 0 I
β
x

 N∑
j=0

q j(x)D
α j y(x)

 (11)
applying equations (3) and (7) to equation (11) gives

y(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(β)

∫ x
0

(x − s)β−1 (12)

×

 N∑
j=0

q j(x)
1

Γ(m j −α j)

∫ s
0

(s − t)m j−α j−1 y(m j )(t)dt

 d s
Substituting equation (4) into equation (12) gives

y(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(β)

∫ x
0

(x − s)β−1 (13)

×

 N∑
j=0

q j(x)
1

Γ(m j −α j)

∫ s
0

(s − t)m j−α j−1
dm j

dtm j
(φ(t)) dtA

 d s

3.1. Method of Solution
Collocating at xi in equation (13) gives

y(xi) = W(xi) +
N∑

j=0

1
Γ(m j −α j)

1
Γ(β)

∫ xi
0

(xi − s)
β−1q j(s)

×

(∫ s
0

(s − t)m j−α j−1
dm j

dtm j
(φ(t)) dt

)
d s A (14)

where

W(xi) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(β)

∫ x
0

(x − s)β−1h(s)d s

2



Ajileye & James / J. Nig. Soc. Phys. Sci. 5 (2023) 1075 3

Simplifying equation (14) gives

φ(xi)A = W(xi)

+


∑N

j=0
1

Γ(m j −α j)
1

Γ(β)

∫ xi
0

(xi − s)β−1q j(s)

×

(∫ s
0

(s − t)m j−α j−1
dm j

dtm j
(φ(t)) dt

)
d s

 A
(15)

Factorizing the values of A from equation (15) gives
φ(xi)−

∑N
j=0

1
Γ(m j −α j)

1
Γ(β)

∫ xi
0

(xi − s)β−1q j(s)

×

(∫ s
0

(s − t)m j−α j−1
dm j

dtm j
(φ(t)) dt

)
d s

 A
= W(xi) (16)

Equation (16) can be in the form

V (xi)A =W(xi) (17)

where

V (xi) = φ(xi)−
N∑

j=0

1
Γ(m j −α j)

1
Γ(β)

∫ xi
0

(xi − s)
β−1q j(s)

(∫ s
0

(s − t)m j−α j−1
dm j

dtm j
(φ(t)) dt

)
d s (18)

and
A = [a0 a1 · · · aN ]T

Multiplying both sides of equation (17) by V−1(xi) gives

A =V−1(xi)W(xi) (19)

Lemma (3.2)
Let y(x) be approximated by (11) and let

L(x) = 0 I
β
x

 N∑
j=0

q j(x)D
α j y(x)

 (20)
If q j(s) = sp j, then

L(x; n) =
Γ(n + 1)Γ(n −α j + p j + 1)

Γ(n −α j + 1)Γ(β + n −α j + p j + 1)

× x
β+n−α j +p j
i A (21)

Proof.
Applying equation (3) and (7) into equation (20) gives

0 I
β
x

 N∑
j=0

q j(x)D
α j y(x)

 = N∑
j=0

1
Γ(m j −α j)

1
Γ(β)

∫ xi
0

(x − s)β−1q j(s)

[∫ s
0

(s − t)m j−α j−1y(m j )(t)dt
]

d s (22)

Substituting (8) into (22) gives

=

N∑
j=0

1
Γ(m j −α j)

1
Γ(β)

∫ x
0

(x − s)β−1S p j (23)

[∫ s
0

(s − t)m j−α j−1
(

Γ(n + 1)
Γ(n − m j + 1)

tn−m j
)

dt
]

d s A

Let s − t = (1 − v)s, then t = vs =⇒
dt
dv

= s =⇒ dt = sdv,
substituting into (23) gives

=

N∑
j=0

Γ(n + 1)
Γ(m j −α j)Γ(n − m j + 1)

1
Γ(β)

∫ x
0

(x − s)β−1S p j (24)

[
S n−α j

∫ 1
0

(1 − v)m j−α j−1V n−m j dt
]

d s A

Simplifying (24), we get

L(x; n) =
Γ(n + 1)Γ(n −α j + p j + 1)

Γ(n −α j + 1)Γ(β + n −α j + p j + 1)
xβ+n−α j +p j A(25)

Lemma (3.3)
Let y(t) be approximated by (9), let

C(x) = 0 I
β
x (h(x)) (26)

if h(s) = sm, then

C(x) =
Γ(m + 1)

Γ(β + m + 1)
xβ+m

Proof.
Applying equation (7) into (26) gives

0 I
β
x (h(x)) =

1
Γ(β)

∫ x
0

(x − s)β−1h(s) d s

Substituting for h(s) gives

=
1

Γ(β)

∫ x
0

(x − s)β−1 smd s

Let x − s = (1 − u)x, s = ux =⇒
d s
du

= x =⇒ d s = xdu.

C(x) =
Γ(m + 1)

Γ(β + m + 1)
xβ+m (27)

Lemma (3.4)
Let y(x) be the solution of (1) and (2) then the numerical result
gives

y(x) = φ(xi)V
−1(xi) W(xi) (28)

where

V (xi) =
Γ(n + 1)Γ(n −α j + p j + 1)

Γ(n −α j + 1)Γ(β + n −α j + p j + 1)
x
β+n−α j +p j
i

and

W(xi) = −
N∑

k=0

y(k)(0)
k!

xki +
Γ(m + 1)

Γ(β + m + 1)
xβ+mi

3



Ajileye & James / J. Nig. Soc. Phys. Sci. 5 (2023) 1075 4

Proof.
Approximate solution of equation (17) is

y(x) = φ(x) A

From equation (19) A =V−1(xi) W(xi)
where

V (xi) =
Γ(n + 1)Γ(n −α j + p j + 1)

Γ(n −α j + 1)Γ(β + n −α j + p j + 1)
x
β+n−α j +p j
i

+
br+n+1Γ(r + 1)

(σ + n + 1) Γ(β + r + 1)
xβ+r

+
Γ(r + σ + n + 2)

(σ + n + 1) Γ(β + r + σ + n + 2)
xβ+r+σ+n+1

Substituting for A in the approximate solution gives the numer-
ical result

y(x) = φ(xi)V
−1(xi) W(xi)

4. Convergence Analysis

In this section, we establish the convergence of the method by
substituting the approximate solution into equation (3.0)

yN (x) = W(x) +
N∑

j=0

1
Γ(m j −α j)

1
Γ(β)

(29)

×

∫ x
0

(x − s)β−1q j(s)
[∫ s

0
(s − t)m j−α j−1y

(m j )
N (t)dt

]
d s

Subtracting (9) from (29) gives

EN (x) = yN (x) − y(x).

Hence

|EN (x)| ≤
1

Γ(β)

∫ x
0

(x − s)β−1
N∑

j=0

1
Γ(m j −α j)

q j(s)∣∣∣∣∣∣
[∫ s

0
(s − t)m j−α j−1 EN (t)dt

]
d s

∣∣∣∣∣∣
Therefore

‖EN (xi)‖∞
‖EN (t)‖∞

≤
1

Γ(β)

∫ xi
0

(x − s)β−1∣∣∣∣∣∣∣∣
 N∑

j=0

1
Γ(m j −α j)

q j(s)
(∫ s

0
(s − t)m j−α j−1dt

)
∣∣∣∣∣∣∣∣ d s

5. Numerical Examples

In this section, we considered two numerical examples to eval-
uate the effectiveness and clarity of the method. A MAPLE
18 program is used to perform the computations. Let yn(x)
and y(x) be the approximate and exact solutions respectively.
ErrorN = |yn(x) − y(x)| .

Example 5.1. [2] Consider multi-order Fractional differential
equation .

D1.5y(x) = −x−1 D0.5y(x) − x0.5y(x) + f (x)

with this condition y
′

(0) = y(0) = 0 and exact solution y(x) =
x3 − x2

f (x) =
[
6x

(
Γ(3.5) + Γ(2.5)
Γ(2.5)Γ(3.5)

+
x2

6

)
− 2

(
Γ(2.5) + Γ(1.5)
Γ(1.5)Γ(2.5)

+
x2

2

)]
x0.5

Solution 1. Comparing with equation (1.1) and Equation (1.2),
β = 1.5,α = 0.5

Using N = 4 for illustration, and applying equation (6) gives

y(x) = W(x) −
1

Γ(1 − 0.5)
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s−1 (30)[∫ s
0

(s − t)1−0.5−1
Γ(n + 1)

Γ(n − 1 + 1)
tn−1dt

]
d s A

−
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s0.5 (sn) d s A

where

W(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 f (s)d s

Substituting (4) into equation (30) gives

φ(x) A = W(x) −
1

Γ(1 − 0.5)
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s−1 (31)[∫ s
0

(s − t)1−0.5−1
Γ(n + 1)

Γ(n − 1 + 1)
tn−1dt

]
d s A

−
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s0.5 (sn) d s A

where

W(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 f (s)d s

Equation (31) can be in the form

τ(x)A =W(x) (32)

where

τ(x) = φ(x) +
1

Γ(1 − 0.5)
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s−1[∫ s
0

(s − t)1−0.5−1
Γ(n + 1)

Γ(n − 1 + 1)
tn−1dt

]
d s

+
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s0.5 (sn) d s

4



Ajileye & James / J. Nig. Soc. Phys. Sci. 5 (2023) 1075 5

collocating at x4 =
[

1
4

2
4

3
4 1

]
and substituting the

initial conditions gives

τ(xi)
∗A = f (xi)∗ (33)

where

τi (x)
∗

=

0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
0.0470157986 0.2239605405 0.1797773130 0.0750069286 0.0274673600
0.1880631945 0.5184447788 0.7543711011 0.6176863534 0.4482932221
0.4231421877 0.9539764130 1.8295669120 2.1838653930 2.3438648990
0.7522527781 1.6010791410 3.5816739880 5.5056804110 7.7368811360

0 1 0 0 0
1 0 0 0 0



f (x)∗ = [0.0000000000 − 0.1047703844 − 0.1366847478 0.3542984804 1.9240064220]

We now solve for the unknown values A making use of matrix
inversion results in equation (33):

y4 =
(

1.365574320288940 × 10−14 − 2.546407529280260 × 10−12 x
−0.999999999275360x2 + 0.999999998413614x3 + 6.644427230639850 × 10−10 x4

)

Example 5.2. [2] Consider multi-order Fractional differential
equation .

D1.5y(x) +
1
x

D0.5y(x) + x
1
2 y(x) = + f (x)

with this condition y
′

(0) = y(0) = 0 and exact solution y(x) =
−x3 + x2

f (x) =
[
2
(
Γ(2.5) + Γ(1.5)
Γ(1.5)Γ(2.5)

+
x2

2

)
− 6x

(
Γ(3.5) + Γ(2.5)
Γ(2.5)Γ(3.5)

+
x2

9

)]
x

1
2

Solution 2. Comparing with equation (1.1) and Equation (1.2),
β = 1.5,α = 0.5

Using N = 4 for illustration,
Applying equation (6) gives

y(x) = W(x) −
1

Γ(1 − 0.5)
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s−1(34)[∫ s
0

(s − t)1−0.5−1
Γ(n + 1)

Γ(n − 1 + 1)
tn−1dt

]
d s A

−
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s0.5 (sn) d s A

Table 1: Exact, approximate and absolute error values for Example 1

x Exact Our methodN=4 errorN=4 error [2]=4
0.0 0.00000000000 1.36557432000e-14 1.3656e-14 5.9232e-13
0.1 -0.900000000e-2 -0.89999999950e-2 5.0000e-12 2.6668e-10
0.2 -0.320000000e-1 -0.31999999980e-1 2.0000e-11 9.6994e-10
0.3 -0.630000000e-1 0.06299999997000 3.0000e-11 1.9604e-09
0.4 -0.960000000e-1 -0.95999999980e-1 2.0000e-11 3.0781e-09
0.5 -0.12500000000 -0.1250000000000 0.0000000 4.1532e-09
0.6 -0.14400000000 -0.1439999999000 1.0000e-10 5.0056e-09
0.7 -0.14700000000 -0.1470000000000 0.0000000 5.4456e-09
0.8 -0.12800000000 -0.1280000001000 1.0000e-10 5.2733e-09
0.9 -0.810000000e-1 -0.81000000160e-1 1.6000e–10 4.2787e-09
1.0 0.00000000000 -2.3555727690e-10 2.3555e-10 2.2421e-09

where

W(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 f (s)d s

Substituting (4) into equation (34) gives

φ(x) A = W(x) −
1

Γ(1 − 0.5)
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s−1(35)[∫ s
0

(s − t)1−0.5−1
Γ(n + 1)

Γ(n − 1 + 1)
tn−1dt

]
d s A

−
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s0.5 (sn) d s A

where

W(x) =
N∑

k=0

y(k)(0)
k!

xk +
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 f (s)d s

Equation (35) can be in the form

τ(x)A =W(x) (36)

where

τ(x) = φ(x) +
1

Γ(1 − 0.5)
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s−1[∫ s
0

(s − t)1−0.5−1
Γ(n + 1)

Γ(n − 1 + 1)
tn−1dt

]
d s

+
1

Γ(1.5)

∫ x
0

(x − s)1.5−1 s0.5 (sn) d s

Collocating at x4 =
[

1
4

2
4

3
4 1

]
and substituting the

initial conditions gives

τ(xi)
∗A = f (xi)∗ (37)

where

τi (x)
∗

=

0.5000000000 2.3817583340 1.9118819450 0.7976779168 0.2921077680
0.7071067812 1.9493225120 2.8363918970 2.3224651180 1.6855566990
0.8660254038 1.9524590850 3.7444893700 4.4696155660 4.7970791030
1.0000000000 2.1283791670 4.7612638900 7.3189233350 10.2849485700

1 0 0 0 0
0‘ 1 0 0 0



f (x)∗ =
[

1.1142040280 0.5139267798 -0.7251261978 -2.5576594460 0 0
]

We now solve for the unknown values A making use of matrix
inversion results in equation (37);

y4 =
(

1.428190898877800 × 10−12 − 2.777014174171200 × 10−10 x
+1.000000001121990x2 − 1.000000001716120x3 + 6.387779194483300 × 10−10 x4

)

5



Ajileye & James / J. Nig. Soc. Phys. Sci. 5 (2023) 1075 6

Table 2: Exact, approximate and absolute error values for Example 2

x Exact Our methodN=4 errorN=4 error [2]=4
0.0 0.00000000000 1.428190899000e-12 1.4281908990000e-12 3.7782e-12
0.1 0.00900000000 0.00899999998200 1.8000000000000e-11 2.4706e-09
0.2 0.03200000000 0.03199999997000 3.0000000000000e-11 1.4306e-08
0.3 0.06300000000 0.06299999997000 3.0000000000000e-11 3.9585e-08
0.4 0.09600000000 0.09599999999000 1.0000000000000e-11 78988e-08
0.5 0.12500000000 0.12499999990000 1.0000000000000e-10 1.2980e-07
0.6 0.14400000000 0.14399999990000 1.0000000000000e-10 1.8590e-07
0.7 0.14700000000 0.14699999980000 2.0000000000000e-10 2.3778e-07
0.8 0.12800000000 0.12799999970000 3.0000000000000e-10 2.7253e-07
0.9 0.08100000000 0.08099999952000 4.8000000000000e-10 2.7386e-07
1.0 0.00000000000 -3.6122208060e-10 3.6122208060000e-10 2.2207e-07

6. Discussion of Results

In this section, we discuss the numerical results obtained by
applying the derived numerical method to the solved examples.
We observed from the result obtained for example 1 as shown
in Table 1 that the approximate solution at N=4 gives y4(x) =
1.365574320288940 × 10−14 − 2.546407529280260 × 10−12 x+
1.000000001121990x2 − 1.000000001716120x3 +
6.387779194483300 × 10−10 x4. The numerical result al-
most converges to the exact solution and produces extremely
small errors. This demonstrated that our method outperformed
the proposed method by Uwaheren et al (2020).

The results of the numerical example 2 in Table 2 shows the
approximate solution at N=4 as y4(x) = 1.428190898877800 ×
10−12 − 2.777014174171200 × 10−10 x + 1.0000001121990x2 −
1.0000001716120x3 + 6.387779194483300 × 10−10 x4. The nu-
merical result converge to the exact solution and give better re-
sult than the method proposed by Uwaheren et al (2020) at the
same value of N. This shows that the numerical method devel-
oped is consistent and converges faster.

7. Conclusion

In this paper, a new numerical method was developed for solv-
ing multi-order fractional differential equations with initial con-
ditions using collocation method. The numerical method de-
rived is consis- tent, efficient and reliable and easy to compute.
Maple code was used to implement the developed method.
Solved numerical examples show that the method is reliable and
suitable for these kind of problems. We also compare our ab-
solute errors with Uwaheren et al. (2020) as shown in Tables 1
and 2. Hence, we safely conclude that our method is preferable
to the existing methods.

References

[1] S.Abbas, & D.Mehdi, “A new operational matrix for solving fractional or-
der differential equations”, Computer and Mathematics with Application
59 (2010) 1326, doi: 10.1016/j.camwa.2009.07.006.

[2] O. A. Uwaheren, A. F. Adebisi & O. A. Taiwo, “Perturbed Collocation
Method For Solving Singular Multi-order Fractional Differential Equa-
tions of Lane-Emden Type”, Journal of the Nigerian Society of Physical
Sciences 3 (2020) 141, https://doi.org/10.46481/jnsps.2020.69.

[3] A. M. Wazwaz & S. M. El-Sayed, “A new modification of the Adomian
decompostion method for linear and nonlinear operators”, App. Math.
Comput. 122 (2001) 393.

[4] R. H. Khan & H. O. Bakodah, “Adomian decomposition method and
its modification for nonlinear Abel’s integral equations”, Computers and
Mathematics with Applications 7 (2013) 2349.

[5] R. C. Mittal & R. Nigam. ”Solution of fractional integro-differential equa-
tions by Adomiandecomposition method”, The International Journal of
Applied Mathematics and Mechanics 2 (2008) 87.

[6] D. A. Gegele, O. P. Evans & D. Akoh, “Numerical solution of higher
order linear Fredholm integro-differential equations”, American Journal
of Engineering Research 3 (2014) 243.

[7] O. A. Agbolade & T. A. Anake, “Solutions of first-order Volterra type
linear integro differential equations by collocation method”, Journal of
Applied Mathematics 4 (2017) 1. doi: 10.1155/2017/1510267.

[8] S. Nemati, P. Lima & Y. Ordokhani, “Numerical method for the mixed
Volterra-Fredholm integral equations using hybrid Legendre function”,
Conference Application of Mathematics 4 (2015) 184.

[9] A. O. Adesanya, Y. A.Yahaya, B. Ahmed & R. O. Onsachi, “Numeri-
cal Solution of Linear integral and Integro-Differential Equations Using
Boubakar Collocation Method”, International Journal of Mathematical
Analysis and Optimization: Theory and Application 2 (2019) 592.

[10] K. Issa & F. Saleh, “Approximate solution of pertubed Volterra Fredholm
integro differential equation by Chebyshev-Galerkin method”, Journal of
Mathematics 17 (2017) 6. doi:10,1155/2017/8213932.

[11] A. H. Bhraway, E. Tohidi & F. Soleymani, “A new Bernoulli matrix
method for solving high order linear and nonlinear Fredholm integro-
differential equations with piecewise interval”, Appl. Math. Comput.
219(2012) 482.

[12] C. Ercan & T. Kharerah, “Solving a class of Volterra integral system by
the differential transform method”, Int. J. Nonlinear Sci. 16 (2013) 87.

[13] M. El-kady & M. Biomy, “Efficient Legendre pseudospectral method for
solving integral and integro differential equation”, Commom Nonlinear
Sci. Numer. Simulat. 15 (2010) 1724.

[14] N. Irfan, S. Kumar & S. Kapoor, “Bernstein Operational Matrix Approach
for Integro-Differential Equation Arising in Control theory ”, Nonlinear
Engineering 3 (2014) 117.

[15] M. K. Shahooth, R. R. Ahmed, U-K. S. Din, W. Swidan, O. K. Al-
Husseini & W. K. Shahooth, “Approximation Solution to Solving Lin-
ear Volterra-Fredholm Integro-Differential Equations of the Second Kind
by Using Bernstein Polynomials Method”, J. Appl. Computat. Math. 5
(2016) 228, DOI:10.4172/2168-9679.1000298.

[16] S. E. Fadugba, “Solution of Fractional Order Equations in the Domain
of the Mellin Transform”, Journal of the Nigerian Society of Physical
Sciences 4 (2019) 138, https://doi.org/10.46481/jnsps.2019.31.

[17] A. Bolandtalat, E. Babolian, H. Jafari, “Numerical solution of multi-order
fractional differential equations by Boubakar polynomial”, Open Phys 14
(2016) 226.

[18] G. Ajileye, A. A. James, A. M. Ayinde & T. Oyedepo, “Collocation Ap-
proach for the Computational Solution Of Fredholm-Volterra Fractional
Order of Integro-Differential Equations”, J. Nig. Soc. Phys. Sci. 4 (2022)
834.

6