EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 4, Article Number 6763
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Employing Cubic B-Spline Functions to Solve Linear
Systems of Volterra Integro-Fractional Differential

Equations with Variable Coefficients
Diar Kh. Abdullah1,∗, Shazad Sh. Ahmed1, Karwan H. F. Jwamer1
1 Mathematics Department, College of Science, University of Sulaimani, Sulaimanyah, Iraq

Abstract. In this study, we present a collocation method based on cubic B-spline functions for
solving systems of Volterra integro-differential equations involving both classical and fractional
derivatives in the Caputo sense (LSVIDEs-CF). The approach begins by dividing the problem
domain into a finite number of subintervals, followed by the construction of cubic B-spline basis
functions within each segment. Control points are introduced as the unknowns in the approxi-
mate numerical solution, which is expressed as a cubic combination of these basis functions. The
given system of VIFDEs-CF is then reduced to a system of algebraic equations, which are ef-
ficiently solved using the Jacobian matrix method. In practice, the integrals are approximated
using the Clenshaw–Curtis quadrature rule. The implementation of this method is supported by
Python software, ensuring effective computational processing. Numerical examples are provided
to demonstrate the efficiency of the proposed method. An itemized version of the algorithm is also
presented to facilitate its implementation.
2020 Mathematics Subject Classifications: 65R20, 65M70, 34A08, 45J05
Key Words and Phrases: System of Volterra integro-fractional differential equations, cubic
B-spline function, collocation method, Jacobian matrix algorithm, Clenshaw-Curtis quadrature
rule

1. Introduction

Integro-differential equations (IDEs) are widely used in mathematical modeling be-
cause they combine integration and differentiation of an unknown function. Scholars such
as Volterra, Fredholm, Malthus, Verhulst, Abel, and Lotka have employed IDEs to study
problems in fields including mathematical biology, physics, and economics [1]. Over the
past several decades, extensive research has been conducted on IDEs, resulting in nu-
merous publications and theoretical advancements. One notable application of IDEs is in

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6763

Email addresses: Diar.Khalid85@gmail.com (D. Kh. Abdullah),
Shazad.ahmed@univsul.edu.iq (Sh. Sh. Ahmed), karwan.jwamer@univsul.edu.iq (K. H. F. Jwamer)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 2 of 40

modeling the dynamic spread of viruses using statistical methods [2, 3]. Fractional integro-
differential equations, in particular, have gained attention due to their ability to represent
memory effects and hereditary behavior in complicated systems, and hence, the search
for correct approximate solutions has developed into one of the central fields of study in
numerical analysis, which led to significant results in the numerical solution of linear types
of these equations [4]. Various kinds of fractional integral equations have been resolved
successfully by building various approaches in recent times [5]. For such complex systems,
reliable numerical methods are necessary to avoid misleading predictions. The flexibility
of the cubic B-spline collocation method is that it can be effectively applied to both linear
and nonlinear FIDEs, arising in viscoelasticity and population control, demonstrating its
usefulness in applications [6]-[7]. Current research has been towards the development of
numerical methods that can provide precise approximations of the solutions of FIDEs.
Among them, the cubic B-spline collocation technique was recently popularized, as cubic
B-splines inherently possess the flexibility and smoothness that would be suitable for ap-
proximating fractional derivative functions [8]. Some researchers have presented extremely
powerful methods to derive efficient and viable approximate solutions of these equations.
Yi et al. presented a graphical illustration of the CAS wavelet method for solving FIDEs
[9]. Quantic B-spline techniques were proposed by Raslan et al. and Fazal et al. see
[10, 11]. Diar et al. illustrated quadratic and linear spline functions, see [7, 12]. While
Mohammad et al. demonstrated the spline collocation method [13]. Abbas et al. in-
troduced an efficient spline technique for solving time-fractional problems [14]. Ghosh et
al. presented an iterative difference scheme [15]. Masti et al. worked on the collocation-
Galerkin method [16]. Khan et al. applied the LADM procedure [17]. Khader et al. used
the Chebyshev pseudo-spectral method [18]. Huang et al. presented the Taylor expansion
method [19]. Yn H et al. demonstrated the Taylor series method [20]. Singh et al. ap-
plied the adaptive mesh method [21], Li et al. explored deep learning techniques [22], and
Shazad worked on the system of linear Volterra integro-fractional differential equations
[23]. In this paper, our focus is on approximating a (SVIDE’s-CF) using a new extended
cubic B-spline based collocation method. The target (SVIDE’s-CF) with the general forms
defined on any closed-bounded interval [a, b] is described by the equations:

Pi(x)Y ′
i(x) + ai0(x)

C
aDσi0

x Yi(x) + ai1(x)Yi(x) = Fi(x) +

m∑
j=0

ωij

∫ x

a
Kij(x, s)

C
aD

βij
s Yj(s) ds.

(1)
Under the following conditions:[

Dki
x Yi(x)

]
x=a

= ϑiki , ∀ki = 0, 1, . . . , µi − 1, and i = 0, 1, . . . ,m. (2)

The variable coefficients Pi(x)(6≡ 0), ai0(x) and ai1(x) ∈ C([a, b],R) and Kij ∈ C(Θ,R),Θ =
{(x, s) : a ≤ s < x ≤ b}, for each i, j = 0, 1, . . . ,m, with fractional orders: σi1 > σi0 > 0
and βim > βi(m−1) > · · · > βi1 > βi0 = 0, and for all i, j = 0, 1, . . . ,m. Furthermore, the
µi = max

{
1,mβ

im

}
for all i = 0, 1, . . . ,m, where mβ

ij−1 < βij ≤ mβ
ij ,m

β
ij = dβije , ωij ∈ R

also for all i, j = 0, 1, . . . ,m. In this work, we examined and assessed a system of



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 3 of 40

Volterra integro-differential equations for classical and fractional orders (SVIDE’s-CF), us-
ing quadrature techniques in conjunction with cubic B-spline functions. The information
was summarized by an algorithm, and we later produced Python software to implement
the method and a few test instances that were resolved with the help of these algorithms.
This paper is structured in the following manner: Section 2 presents background informa-
tion on fractional calculus and highlights its key features. Section 3 discusses the materials
used and introduces cubic B-spline functions along with new lemmas. Section 4 explores
the derivation of the proposed method. Section 5 provides the algorithms. Section 6 com-
pares the numerical results with those from existing methods. Finally, Section 7 offers a
concise summary of the study’s findings.

2. Basic definitions of fractional derivatives

This section discusses fractional calculus and its properties and defines the numerical
integration as the Clenshaw-Curtis quadrature rule:

Definition 1. ([3], [24]). The fractional operators aJ
α
x and R

a Dα
x of order α for a func-

tion Y(x) are called the Riemann-Liouville fractional integral and differential operators,
respectively, for all n− 1 < α ≤ n (n ∈ Z+, α ∈ R+).

aJ
α
x Y(x) =

1

Γ(α)

∫ x

a

Y(s)

(x− s)1−α
ds; aJ

0
xY(x) = Y(x), x > a (3)

R
a Dα

xY(x) =
dn

dxn
(
aJ

n−α
x Y(x)

)
=

1

Γ(n− α)

dn

dxn

∫ x

a

Y(s)

(x− s)α+1−n
ds, x > a (4)

Definition 2. ([3], [24]). The fractional operator C
aDα

x of order α for a function Y(x) is
called the Caputo fractional differential operator, for all n−1 < α ≤ n (n ∈ Z+, α ∈ R+),
and is defined as:

C
aDα

xY(x) = aJ
n−α
x

(
dn

dxn
Y(x)

)
=

1

Γ(n− α)

∫ x

a

Y(n)(s)

(x− s)α+1−n
ds, x > a (5)

The following are some interesting properties of the operator C
aDα

x :

• For n = α, we have
C
aDα

xY(x) =
dn

dxn
Y(x).

• The operator C
aD

α
x is linear, i.e., for any ϕ1, ϕ2 ∈ R,

C
aDα

x [ϕ1Y1(x)± ϕ2Y2(x)] = ϕ1
C
aDα

t Y1(x)± ϕ2
C
aDα

xY2(x).

• For any constant ϕ ∈ R,
C
aDα

xϕ = 0.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 4 of 40

• If Y(x) = (x− a)p, then:

C
aDα

x (x− a)p =


Γ(p+ 1)

Γ(p− α+ 1)
(x− a)p−α, if p ∈ N, p ≥ dαe or p /∈ N, p > dαe − 1

0, if p ∈ {0, 1, 2, . . . , dαe − 1}

Definition 3. ([12], [25], [26]). Clenshaw and Curtis established a method for approxi-
mating a definite integral by expressing the integrand as a finite Chebyshev series. This
approach is especially effective when applied to integral equations.

∫ 1

−1
Y(x) dx ≈

N∑
r=0, even

′′

(
2

N

N∑
k=0

cos

(
rkπ

N

)
Y
(
cos

(
kπ

N

)))
, k = 0, 1, . . . ,N . (6)

Remark 2.1 ([26]). The following observations can be made:

• The notation
∑′′ indicates that the first and last terms in the summation are to be

halved.

• The transformation
x =

b− a

2
t+

b+ a

2

maps x ∈ [−1, 1] to t ∈ [a, b].

3. Materials and cubic B-spline function

In this section, we will define the cubic B-spline function, the B-spline curve and knot
vectors, starting with the definition of a B-spline of degree k.

Definition 4. ([27], [28]). Let TN = {x0, x1, . . . , xN} be a uniform or non-uniform
partition of the interval [a, b]. The k-degree B-spline basis function Bk

r (x), for r ≥ 0, is
defined recursively as:

Bk
r (x) =

x− xr
xr+k − xr

Bk−1
r (x) +

xr+k+1 − x

xr+k+1 − xr+1
Bk−1
r+1 (x) (7)

Remark 1 ([29]). We note the following facts about the B-spline basis functions:

• Local support: Bk
r (x) = 0 for all x /∈ [xr, xr+k+1).

• Nonnegativity: Bk
r (x) > 0 for all x ∈ [xr, xr+k+1).

Remark 2 ([27], [30]). In the context of B-spline curves, the knot vector T =
{x0, x1, x2, . . . , xm} plays a crucial role in determining the shape of the curve.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 5 of 40

Definition 5. ([27], [30]). An open uniform knot vector is a knot vector that contains k
equal knot values at each end. The parametric knot values xr for an open curve are given
by:

xr =


0, r < k + 1

r − k, k + 1 ≤ r ≤ n

n− k + 1, r > n

(8)

where 0 ≤ r ≤ n+ k + 1, and the parameter x lies in the range 0 ≤ x ≤ n− k + 1.

Definition 6. ([27]). A B-spline curve is a piecewise-defined polynomial curve. Mathe-
matically, a B-spline curve B(k,n)(x) of degree k, with control points {ϑ0, ϑ1, . . . , ϑn} and
a non-decreasing knot vector T = {x0, x1, x2, . . . , xm}, is defined as:

B{k,n}(x) =
n∑

z=0

ϑz Bk
z (x), x0 ≤ x ≤ xm (9)

Here, ϑz are the control points that influence the shape of the curve, and Bk
z (x) are

the B-spline basis functions of degree k, which are defined recursively and determine the
influence of each control point on the curve for any parameter value x.

Lemma 1. General expression of a cubic B-spline curve. Let ϑ0, ϑ1, ϑ2, and ϑ3 be four
consecutive control points of a third-degree B-spline curve defined on the interval [a, b].
Then, the curve is given by:

B{k,n}(x?) = ϑ0

(
b− x?

b− a

)3

+ 3ϑ1

(
x? − a

b− a

)(
b− x?

b− a

)2

+ 3ϑ2

(
x? − a

b− a

)2(b− x?

b− a

)
+ ϑ3

(
x? − a

b− a

)3

, x? ∈ [a, b] (10)

Proof. First, we convert the given knots xr ∈ T , defined on the original interval 0 ≤ x ≤
n−k+1 as in (8), to a new knot vector x?r ∈ T ? defined on a new interval x? ∈ [a, b]. This
is achieved by applying the following linear transformation to each element of the original
knot vector T = {x0, x1, x2, . . . , xm}, yielding a new knot vector T ? = {x?0, x?1, . . . , x?m}:

x?r = a+

(
b− a

n− k + 1

)
xr (11)

By applying the transformation defined in (11) to each case in (8), we obtain:

x?r =


a, r < k + 1

a+

(
b− a

n− k + 1

)
(r − k), k + 1 ≤ r ≤ n

b, r > n

(12)

Where:



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 6 of 40

• For r < k + 1, xr = 0 implies x?r = a

• For k + 1 ≤ r ≤ n, xr = r − k implies x?r = a+

(
b− a

n− k + 1

)
(r − k)

• For r > n, xr = n− k + 1 implies x?r = b

Substituting k = 3 (cubic B-spline) into the recurrence relation from (7), we obtain the
cubic B-spline basis functions used in the expression above.

B3
r(x

?) =
x? − x?r
x?r+3 − x?r

B2
r(x

?) +
x?r+4 − x?

x?r+4 − x?r+1

B2
r+1(x

?) (13)

We substitute the expressions for B2
r(x

?) and B2
r+1(x

?), which are computed using the
first and zero degree B-spline basis functions. After simplification and substitution, we
obtain:

B3
r(x

?) =

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

(x? − x?r)
3

(x?r+3 − x?r)(x
?
r+2 − x?r)(x

?
r+1 − x?r)

, if x?r ≤ x? < x?r+1

(x? − x?r)
2(x?r+2 − x?)

(x?r+3 − x?r)(x
?
r+2 − x?r)(x

?
r+2 − x?r+1)

+
(x? − x?r)(x

?
r+3 − x?)(x? − x?r+1)

(x?r+3 − x?r)(x
?
r+3 − x?r+1)(x

?
r+2 − x?r+1)

+
(x?r+4 − x?)(x? − x?r+1)

2

(x?r+4 − x?r+1)(x
?
r+3 − x?r+1)(x

?
r+2 − x?r+1)

if x?r+1 ≤ x? < x?r+2

(x?r+3 − x?)2(x? − x?r)

(x?r+3 − x?r)(x
?
r+3 − x?r+1)(x

?
r+3 − x?r+2)

+
(x?r+4 − x?)(x? − x?r+1)(x

?
r+3 − x?)

(x?r+4 − x?r+1)(x
?
r+3 − x?r+1)(x

?
r+3 − x?r+2)

+
(x?r+4 − x?)2(x? − x?r+2)

(x?r+4 − x?r+1)(x
?
r+4 − x?r+2)(x

?
r+3 − x?r+2)

if x?r+2 ≤ x? < x?r+3

(x?r+4 − x?)3

(x?r+4 − x?r+1)(x
?
r+4 − x?r+2)(x

?
r+4 − x?r+3)

, if x?r+3 ≤ x? < x?r+4

0, otherwise

(14)

Now, if we set k = 3 and n = 3, then from (12), we obtain the open uniform knot
vector: T ? = {a, a, a, a, b, b, b, b}. The parameter x? lies in the range: x? ∈ [a, b].

Since the knot vector is uniform and open, x?r ∈ T ? for r = 0, . . . , n. Substituting into
(9) and using the basis functions from (13), we get:

B{k,n}(x?) =
3∑

z=0

ϑz B3
z(x

?)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 7 of 40

= ϑ0B3
0(x

?) + ϑ1B3
1(x

?) + ϑ2B3
2(x

?) + ϑ3B3
3(x

?) (15)

That is:
B{k,n}(x?) =

ϑ0 ·

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

(x? − x?0)
3

(x?3 − x?0)(x
?
2 − x?0)(x

?
1 − x?0)

, if x?0 ≤ x? < x?1

(x? − x?0)
2(x?2 − x?)

(x?3 − x?0)(x
?
2 − x?0)(x

?
2 − x?1)

+
(x? − x?0)(x

?
3 − x?)(x? − x?1)

(x?3 − x?0)(x
?
3 − x?1)(x

?
2 − x?1)

+
(x?4 − x?)(x? − x?1)

2

(x?4 − x?1)(x
?
3 − x?1)(x

?
2 − x?1)

if x?1 ≤ x? < x?2

(x?3 − x?)2(x? − x?0)

(x?3 − x?0)(x
?
3 − x?1)(x

?
3 − x?2)

+
(x?4 − x?)(x? − x?1)(x

?
3 − x?)

(x?3 − x?1)
2(x?3 − x?2)

+
(x?4 − x?)2(x? − x?2)

(x?4 − x?1)(x
?
4 − x?2)(x

?
3 − x?2)

if x?2 ≤ x? < x?3

(x?4 − x?)3

(x?4 − x?1)(x
?
4 − x?2)(x

?
4 − x?3)

, if x?3 ≤ x? < x?4

0, otherwise

+ϑ1 ·

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

(x? − x?1)
3

(x?4 − x?1)(x
?
3 − x?1)(x

?
2 − x?1)

, if x?1 ≤ x? < x?2

(x? − x?1)
2(x?3 − x?)

(x?4 − x?1)(x
?
3 − x?1)(x

?
3 − x?2)

+
(x? − x?1)(x

?
4 − x?)(x? − x?2)

(x?4 − x?1)(x
?
4 − x?2)(x

?
3 − x?2)

+
(x?5 − x?)(x? − x?2)

2

(x?5 − x?2)(x
?
4 − x?2)(x

?
3 − x?2)

if x?2 ≤ x? < x?3

(x?4 − x?)2(x? − x?1)

(x?4 − x?1)(x
?
4 − x?2)(x

?
4 − x?3)

+
(x?5 − x?)(x? − x?2)(x

?
4 − x?)

(x?5 − x?2)(x
?
4 − x?2)(x

?
4 − x?3)

+
(x?5 − x?)2(x? − x?3)

(x?5 − x?2)(x
?
5 − x?3)(x

?
4 − x?3)

if x?3 ≤ x? < x?4

(x?5 − x?)3

(x?5 − x?2)(x
?
5 − x?3)(x

?
5 − x?4)

, if x?4 ≤ x? < x?5

0, otherwise



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 8 of 40

+ϑ2 ·

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

(x? − x?2)
3

(x?5 − x?2)(x
?
4 − x?2)(x

?
3 − x?2)

, if x?2 ≤ x? < x?3

(x? − x?2)
2(x?4 − x?)

(x?5 − x?2)(x
?
4 − x?2)(x

?
4 − x?3)

+
(x? − x?2)(x

?
5 − x?)(x? − x?3)

(x?5 − x?2)(x
?
5 − x?3)(x

?
4 − x?3)

+
(x?6 − x?)(x? − x?3)

2

(x?6 − x?3)(x
?
5 − x?3)(x

?
4 − x?3)

if x?3 ≤ x? < x?4

(x?5 − x?)2(x? − x?2)

(x?5 − x?2)(x
?
5 − x?3)(x

?
5 − x?4)

+
(x?6 − x?)(x? − x?3)(x

?
5 − x?)

(x?6 − x?3)(x
?
5 − x?3)(x

?
5 − x?4)

+
(x?6 − x?)2(x? − x?4)

(x?6 − x?3)(x
?
6 − x?4)(x

?
5 − x?4)

if x?4 ≤ x? < x?5

(x?6 − x?)3

(x?6 − x?3)(x
?
6 − x?4)(x

?
6 − x?5)

, if x?5 ≤ x? < x?6

0, otherwise

+ϑ3 ·

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

(x? − x?3)
3

(x?6 − x?3)(x
?
5 − x?3)(x

?
4 − x?3)

, if x?3 ≤ x? < x?4

(x? − x?3)
2(x?5 − x?)

(x?6 − x?3)(x
?
5 − x?3)(x

?
5 − x?4)

+
(x? − x?3)(x

?
6 − x?)(x? − x?4)

(x?6 − x?3)(x
?
6 − x?4)(x

?
5 − x?4)

+
(x?7 − x?)(x? − x?4)

2

(x?7 − x?4)(x
?
6 − x?4)(x

?
5 − x?4)

if x?4 ≤ x? < x?5

(x?6 − x?)2(x? − x?3)

(x?6 − x?3)(x
?
6 − x?4)(x

?
6 − x?5)

+
(x?7 − x?)(x? − x?4)(x

?
6 − x?)

(x?7 − x?4)(x
?
6 − x?4)(x

?
6 − x?5)

+
(x?7 − x?)2(x? − x?5)

(x?7 − x?4)(x
?
7 − x?5)(x

?
6 − x?5)

if x?5 ≤ x? < x?6

(x?7 − x?)3

(x?7 − x?4)(x
?
7 − x?5)(x

?
7 − x?6)

, if x?6 ≤ x? < x?7

0, otherwise

Through the utilization of T ? and carrying out subsequent computations, we obtain (10):

B{C,3}(x?) =

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

ϑ0(b− x?)3

(b− a)3
+

3ϑ1(b− x?)2(x? − a)

(b− a)3

+
3ϑ2(x

? − a)2(b− x?)

(b− a)3
+

ϑ3(x
? − a)3

(b− a)3

if a ≤ x? ≤ b

0, otherwise

(16)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 9 of 40

The proof is done.

Lemma 2. Let B{C,3}(x?) represent a cubic B-spline curve defined over the closed bounded
interval [a, b] with control points ϑ0, ϑ1, ϑ2, and ϑ3 being four consecutive control points.
The α-Caputo fractional derivative of B{C,3}(x?), for α ∈ (0, 1], is expressed as:

C
aDα

x?B{C,3}(x?) =
3(x? − a)1−α

(Nh)3 Γ(4− α)

[
ϑ0

(
2(Nh)(3− α)(x? − a)− (Nh)2(2− α)(3− α)

− 2(x? − a)2
)

+ ϑ1

(
(Nh)2(2− α)(3− α)− 4(Nh)(3− α)(x? − a)

+ 6(x? − a)2
)

+ ϑ2

(
(Nh)(3− α)(x? − a)− 6(x? − a)2

)
+ 2ϑ3(x

? − a)2
]

(17)
Proof: Use a general expression of a cubic B-spline curve defined on the interval

x? ∈ [a, b]. From Equation (16), it is given by:

B{C,3}(x?) = ϑ0

(
b− x?

b− a

)3

+ 3ϑ1

(
x? − a

b− a

)(
b− x?

b− a

)2

+ 3ϑ2

(
x? − a

b− a

)2(b− x?

b− a

)
+ ϑ3

(
x? − a

b− a

)3

(18)

We simplify the given expression as follows:

B{C,3}(x?) =
ϑ0

(b− a)3
(b− x?)3 +

3ϑ1

(b− a)3
(x? − a)(b− x?)2

+
3ϑ2

(b− a)3
(x? − a)2(b− x?) +

ϑ3

(b− a)3
(x? − a)3

(19)

From the properties of the Caputo fractional derivative in Definition 2, for α ∈ (0, 1],
we can express the cubic B-spline curve for all x? ∈ [a, b] as:



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 10 of 40

C
aDα

x?B{C,3}(x?) =
ϑ0

(b− a)3

[
−3(b− a)

Γ(2− α)
(x? − a)1−α +

6(b− a)

Γ(3− α)
(x? − a)2−α

− 6

Γ(4− α)
(x? − a)3−α

]

+
3ϑ1

(b− a)3

[
(b− a)2

Γ(2− α)
(x? − a)1−α − 4(b− a)

Γ(3− α)
(x? − a)2−α

+
6

Γ(4− α)
(x? − a)3−α

]

+
3ϑ2

(b− a)3

[
2(b− a)

Γ(3− α)
(x? − a)2−α − 6

Γ(4− α)
(x? − a)3−α

]

+
6ϑ3

(b− a)3Γ(4− α)
(x? − a)3−α

(20)

A cubic B-spline curve defined over the closed, bounded interval [a, b] with control points
ϑ0, ϑ1, ϑ2, and ϑ3 - where h is the step size and N is the number of iterations - therefore,
(17) yields the desired fractional derivative expression.

4. Derivation of the method

The main aim of this section is to apply the cubic B-spline collocation method in order
to introduce a numerical approximation technique to solve the problem (1), under the
initial condition given by (2). The method expresses the solution as a cubic combination
of B-spline basis functions.

Let a = x?0 < x?1 < x?2 < · · · < x?N−1 < x?N = b, be a uniform mesh partition of the
domain [a, b] defined by the knots x?r , with

h =
x?r+1 − x?r

N
=

b− a

N
, r = 0, 1, . . . ,N − 1.

For all x? ∈ [a, b], let Yi(x) ≈ B{C,3}
i (x?) for each i = 0, 1, . . . ,m, using the collocation

technique with cubic B-spline basis functions to find an approximate solution B{C,3}
i (x?),

given by:

B{C,3}
i (x?) =

n∑
z=0

ϑiz Bk
iz(x

?), i = 0, 1, . . . ,m (21)

In order to solve the approximate solution of (21), we determine the control points ϑiz

for all i = 0, 1, . . . ,m and z = 0, . . . , n.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 11 of 40

• ϑi0 is determined by the initial condition (2).

• ϑi1 is determined from

ϑi1 =
Nh

3

dB{C,3}
i1 (x?)

dx?

∣∣∣∣∣
x?=a

+ ϑi0,

where

dB{C,3}
i1 (x?)

dx?

∣∣∣∣∣
x?=a

≈ Y ′
i(x

?)
∣∣
x?=a

=
1

Pi(x?)

[
Fi(x

?) +

m∑
j=0

ωij

∫ x?

a
Kij(x

?, s) Ca D
βij
s Yj(s) ds

− ai0(x
?) Ca D

σi0
x? Yi(x

?)− ai1(x
?)Yi(x

?)

]∣∣∣∣∣
x?=a

• ϑi2 is determined from

ϑi2 =
(Nh)2

6

d2B{C,3}
i1 (x?)

dx?2

∣∣∣∣∣
x?=a

− ϑi0 + 2ϑi1,

where the second derivative is approximated using a forward difference formula:

d2B{C,3}
i1 (x?)

dx?2

∣∣∣∣∣
x?=a

≈ Y ′′
i (x

?
r) ≈

Y(p)
i (x?r+2)− 2Y(p)

i (x?r+1) + Yi(x
?
r)

h2

In this case, Yi(x
?
r) is known from the initial condition (2), and Y(p)

i (x?r+2) and
Y(p)
i (x?r+1) are determined from the first-order B-spline [31].

• ϑi3 is determined from the system of (VIDE’s-CF) in (1).

Let the unknown function B{C,3}
i (x?) be approximated using cubic B-spline interpolation

of degree 3, as given in (21). Substituting this approximation into (1), we obtain:

Pi(x
?)

d

dx?

(
n∑

z=0

ϑiz Bk
iz(x

?)

)
+ ai0(x

?) CaD
σi0
x?

(
n∑

z=0

ϑiz Bk
iz(x

?)

)

+ ai1(x
?)

n∑
z=0

ϑiz Bk
iz(x

?)

= Fi(x
?) +

m∑
j=0

ωij

∫ x?

a
Kij(x

?, s) CaD
βij
s

(
n∑

z=0

ϑjz Bk
jz(s)

)
ds

(22)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 12 of 40

Consequently, applying the linearity property of the Caputo fractional derivative from
Definition 2, along with the standard derivative, we derive the following expression:

Pi(x
?)

n∑
z=0

ϑiz
d

dx?

[
Bk
iz(x

?)
]
+ ai0(x

?)
n∑

z=0

ϑiz
C
aD

σi0
x?

[
Bk
iz(x

?)
]

+ ai1(x
?)

n∑
z=0

ϑiz Bk
iz(x

?)

= Fi(x
?) +

m∑
j=0

ωij

∫ x?

a
Kij(x

?, s)
n∑

z=0

ϑjz
C
aD

βij
s

[
Bk
jz(s)

]
ds

(23)
By defining x? = x?r+1 and analyzing the collocation points, we employ a cubic B-spline
function (k = 3, n = 3) within the interval [x?r , x

?
r+1]. We then apply the fractional

derivative in the Caputo sense. From (16), the fractional orders σi0, βij ∈ (0, 1] for all
i, j = 0, 1, . . . ,m, yielding results for each r = 0, 1, . . . ,N − 1.

Pi(x
?
r+1)

{
−3ϑi0(b− x?r+1)

2

(Nh)3
+ 3ϑi1

[−2(x?r+1 − a)(b− x?r+1)

(Nh)3
+

(b− x?r+1)
2

(Nh)3

]

+3ϑi2

[−(x?r+1 − a)2

(Nh)3
+

2(b− x?r+1)(x
?
r+1 − a)

(Nh)3

]
+ 3ϑi3 ·

(x?r+1 − a)2

(Nh)3

}

+ai0(x
?
r+1) ·

(x?r+1 − a)1−σi0

(Nh)3 Γ(4− σi0)
·

{
ϑi0

[
− 3(Nh)2(2− σi0)(3− σi0)

+6(Nh)(3− σi0)− 6(x?r+1 − a)2
]
+ ϑi1

[
(Nh)2(2− σi0)(3− σi0)

−4(Nh)(3− σi0)(x
?
r+1 − a) + 6(x?r+1 − a)2

]
+3ϑi2

[
2(Nh)(3− σi0)(x

?
r+1 − a)− 6(x?r+1 − a)2

]
+ 6ϑi3(x

?
r+1 − a)2

}

+ai1(x
?
r+1) ·

{
ϑi0

(
b− x?r+1

Nh

)3

+ 3ϑi1

(
x?r+1 − a

Nh

)(
b− x?r+1

Nh

)2

+3ϑi2

(
x?r+1 − a

Nh

)2(b− x?r+1

Nh

)
+ ϑi3

(
x?r+1 − a

Nh

)3
}

= Fi(x
?
r+1) +

m∑
j=1

ωij

{
r−1∑
q=0

∫ x?
q+1

x?
q

Kij(x
?
r+1, s) ·

(s− a)1−βij

(Nh)3 Γ(4− βij)
·

[
ϑj0

(
−3(Nh)2(2− βij)(3− βij) + 6(Nh)(3− βij)− 6(s− a)2

)
+3ϑj1

(
(Nh)2(2− βij)(3− βij)− 4(Nh)(3− βij)(s− a) + 6(s− a)2

)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 13 of 40

+3ϑj2

(
2(Nh)(3− βij)(s− a)− 6(s− a)2

)
+ 6ϑj3(s− a)2

]
ds

}

∫ x?
r+1

x?
r

Kij(x
?
r+1, s) · (s− a)1−βij

(Nh)3 Γ(4− βij)
·

[
ϑj0

(
−3(Nh)2(2− βij)(3− βij) + 6(Nh)(3− βij)− 6(s− a)2

)
+3ϑj1

(
(Nh)2(2− βij)(3− βij)− 2(s− a)− 4(Nh)(3− βij)(s− a)

)
+6(s− a)2 + 3ϑj2

(
2(Nh)(3− βij)(s− a)− 6(s− a)2

)
+ 6ϑj3(s− a)2

]
ds (24)

The cubic B-spline function over the interval [x?r , x
?
r+1] is optimized to simplify its

representation and promote efficient solution methods in practice, these integrals must
be approximated using the Clenshaw-Curtis quadrature rule [12, 25, 26], hence, (25) is
applicable for r = 0, 1, . . . ,N − 1 and i, j = 0, 1, . . . ,m.

Vσi0
i,x?

r+1
· ϑi3 +Wσi0

i,x?
r+1

= Fi(x
?
r+1) +

m∑
j=1

ωij

{
Er
ij

}
+

m∑
j=1

ωij

{
Ěr
ij

}
· ϑj3 + Gr

ij + Ǧr
i · ϑ03 (25)

To obtain the unique solution of the system described in (26), we require three additional
equations. These are obtained by computing the coefficients ϑi0, ϑi1, and ϑi2. Once these
are determined, the remaining coefficient ϑi3 can be calculated by solving the following
linear system of algebraic equations of size (m+ 1)× (m+ 1):

A ·B = C (26)

A =



Vσ00
0,x?

r+1
− Ǧr

0 −ω01Ěr
01 −ω02Ěr

02 · · · −ω0mĚr
0m

−Ǧr
1 Vσ10

1,x?
r+1

− ω11Ěr
11 −ω12Ěr

12 · · · −ω1mĚr
1m

−Ǧr
2 −ω21Ěr

21 Vσ20
2,x?

r+1
− ω22Ěr

22 · · · −ω2mĚr
2m

...
...

... . . . ...
−Ǧr

m −ωm1Ěr
m1 −ωm2Ěr

m2 · · · Vσm0
m,x?

r+1
− ωmmĚr

mm


(27)

B =


ϑ03

ϑ13

ϑ23
...

ϑm3

 (28)

C =


L0
r

L1
r

L2
r
...

Lm
r

 (29)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 14 of 40

where,

• Vσi0
i,x?

r+1
=

3Pi(x
?
r+1)(x

?
r+1 − a)2

(Nh)3
+

6 ai0(x
?
r+1)(x

?
r+1 − a)3−σi0

(Nh)3Γ(4− σi0)
+ ai1(x

?
r+1) ·

(
x?r+1 − a

Nh

)3

.

• Wσi0
i,x?

r+1
= Pi(x

?
r+1)

{
−3ϑi0(b− x?r+1)

(Nh)3
+ 3ϑi1

[−2(x?r+1 − a)(b− x?r+1)

(Nh)3
+

(b− x?r+1)
2

(Nh)3

]

+3ϑi2

[−(x?r+1 − a)2

(Nh)3
+

2(b− x?r+1)(x
?
r+1 − a)

(Nh)3

]}

+ai0(x
?
r+1) ·

(x?r+1 − a)1−σi0

(Nh)3Γ(4− σi0)
·

[
ϑi0

(
−3(Nh)2(2− σi0)

)
(3− σi0) + 6(Nh)(3− σi0)− 6(x?r+1 − a)2 + ϑi1

(
(Nh)2

)
(2− σi0)(3− σi0)− 4(Nh)(3− σi0)(x

?
r+1 − a) + 6(x?r+1 − a)2

+3ϑi2

(
2(Nh)(3− σi0)(x

?
r+1 − a)− 6(x?r+1 − a)2

) ]

+ai1(x
?
r+1) ·

[
ϑi0

(
b− x?r+1

Nh

)3

+ 3ϑi1

(
x?r+1 − a

Nh

)(
b− x?r+1

Nh

)2

+3ϑi2

(
x?r+1 − a

Nh

)2(b− x?r+1

Nh

)]
.

• Er
ij =

r−1∑
q=0

N∑
ξ=0

λξ Kij(x
?
r+1, sξ) ·

(x?q+1 − x?q)(sξ − a)1−βij

2(Nh)3Γ(4− βij)
·


ϑj0

[
−3(Nh)2(2− βij)(3− βij) + 6(Nh)(3− βij)− 6(sξ − a)2

]
+ 3ϑj1

[
(Nh)2(2− βij)(3− βij)− 4(Nh)(3− βij)(sξ − a)− 2(sξ − a) + 6(sξ − a)2

]
+ 3ϑj2

[
2(Nh)(3− βij)(sξ − a)− 6(sξ − a)2

]


+
N∑
ξ=0

λξ Kij(x
?
r+1, sξ) ·

(x?r+1 − x?r)(sξ − a)1−βij

2(Nh)3Γ(4− βij)
·


ϑj0

[
−3(Nh)2(2− βij)(3− βij) + 6(Nh)(3− βij)− 6(sξ − a)2

]
+ 3ϑj1

[
(Nh)2(2− βij)(3− βij)− 4(Nh)(3− βij)(sξ − a)− 2(sξ − a) + 6(sξ − a)2

]
+ 3ϑj2

[
2(Nh)(3− βij)(sξ − a)− 6(sξ − a)2

]
 .



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 15 of 40

• Ěr
ij =

r−1∑
q=0

N∑
ξ=0

3λξ Kij(x
?
r+1, sξ) (x

?
q+1 − x?q) (sξ − a)3−βij

(Nh)3 Γ(4− βij)

+

N∑
ξ=0

3λξ Kij(x
?
r+1, sξ) (x

?
r+1 − x?r) (sξ − a)3−βij

(Nh)3 Γ(4− βij)

• Gr
i =

r−1∑
q=0

N∑
ξ=0

λξ Ki0(x
?
r+1, sξ)·

(x?q+1 − x?q)

2
·

[
ϑ00

(
b− sξ
Nh

)3

+3ϑ01

(
sξ − a

Nh

)(
b− sξ
Nh

)2

+ 3ϑ02

(
sξ − a

Nh

)2(b− sξ
Nh

)]
+

N∑
ξ=0

λξ Ki0(x
?
r+1, sξ) ·

(x?r+1 − x?r)

2
·

[
ϑ00

(
b− sξ
Nh

)3

+ 3ϑ01

(
sξ − a

Nh

)(
b− sξ
Nh

)2

+ 3ϑ02

(
sξ − a

Nh

)2(b− sξ
Nh

)]

• Ǧr
i =

r−1∑
q=0

N∑
ξ=0

λξ Ki0(x
?
r+1, sξ) ·

(x?q+1 − x?q)

2

(
sξ − a

Nh

)3

+

N∑
ξ=0

λξ Ki0(x
?
r+1, sξ) ·

(x?r+1 − x?r)

2

(
sξ − a

Nh

)3

Li
r = Fi(x

?
r+1) +

m∑
j=1

ωij Er
ij −Wσi0

i,x?
r+1

+ Gr
ij (30)

Where,

t?ξ =

N ′′∑
ξ=0

cos

(
ξπ

N

)
(31)

are the nodes evaluated to correspond to the extrema of the Chebyshev polynomials on
the interval [−1, 1].
The mapped node sξ is calculated by:

sξ =
x?q+1 − x?q

2
t?ξ +

x?q+1 + x?q
2

(32)

The weights λξ are coefficients that multiply the function values at the nodes to approxi-
mate the integral, given by:

λξ =
2

N

N ′′∑
r=0
even

cos

(
rξπ

N

)
(33)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 16 of 40

5. Algorithms

In this section, we present a step-by-step algorithmic description of the proposed ap-
proach. To facilitate practical implementation, an itemized version of the algorithm is
provided below.

Main Algorithm: LSVIFDE’s (Using Cubic B-spline)

Input Data:

• Domain bounds a, b, and total number of steps N .

• Number of equations: m+ 1.

• Functions: Pi(x
?), ai0(x?), ai1(x?), Fi(x

?).

• Parameters: ωij , Kij(x
?, s), σi0, βij , for all i, j = 0, 1, . . . ,m.

Output:

Vector B from (28), containing control points ϑi3, i = 0, 1, . . . ,m.

Steps:

(i) Construct vectors B and C of size m+ 1, and matrix A of size (m+ 1)× (m+ 1).

(ii) Compute the step size:

h =
b− a

N
, with N ∈ N

and define the partition points:

x?r+1 = a+ (r + 1)h, for r = 0, 1, . . . ,N − 1

(iii) Approximate:
dB{C,3}

i (x?)

dx?

∣∣∣∣∣
x?=a

≈ Y ′
i(a)

using the expression:

1

Pi(a)

Fi(a) +

m∑
j=0

ωij

∫ a

a
Kij(a, s)

C
aD

βij
s Yj(s) ds− ai0(a)

C
aD

σi0
a Yi(a)− ai1(a)Yi(a)


(iv) Compute:

ϑi0 = from given initial condition (2)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 17 of 40

ϑi1 =
Nh

3
·
dB{C,3}

i (x?)

dx?

∣∣∣∣∣
x?=a

+ ϑi0

ϑi2 =
(Nh)2

6
·
d2B{C,3}

i (x?)

d(x?)2

∣∣∣∣∣
x?=a

− ϑi0 + 2ϑi1

(v) Use forward difference to approximate:

d2B{C,3}
i (x?)

d(x?)2

∣∣∣∣∣
x?=a

≈
Y(p)
i (x?r+2)− 2Y(p)

i (x?r+1) + Yi(x
?
r)

h2

(vi) Evaluate Yi(x
?
r) from (2), and Y(p)

i (x?r+1), Y
(p)
i (x?r+2) using first-order B-spline.

(vii) Compute entries of C using (29):

C =


L0
r

L1
r
...

Lm
r


(viii) Construct the matrix A and modify its entries using the computed ϑi0, ϑi1, and ϑi2.

(ix) Solve the system A ·B = C, using Clenshaw-Curtis quadrature for numerical inte-
gration.

Output:

B =
[
ϑ03 ϑ13 ϑ23 . . . ϑm3

]T
First Procedure: Normalize Finding Cubic B-Spline for Each Interval
(NFCBI)

We first perform all steps outlined in the Main Algorithm (LSVIFDE’s), and then
proceed with the following additional steps:

(i) For r = 0, 1, . . . ,N − 1, set:

x? = x?r+1 = a+ (r + 1)h, n = 3, k = 3

(ii) For each i = 0, 1, . . . ,m, evaluate the approximate cubic B-spline solution as:

B{C,3}
i (x?r+1) ≈

n∑
z=0

ϑiz Bk
iz(x

?
r+1)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 18 of 40

Output: The set of approximate solutions for each function:

B{C,3}
0 (x?r+1), B{C,3}

1 (x?r+1), . . . , B{C,3}
m (x?r+1)

represent the cubic B-spline approximation to each Yi(x) over the interval [x?r , x?r+1].

II. Second Procedure: First Control Point Cubic B-Spline (FCPCBS)

We begin by performing all steps from the Main Algorithm (LSVIFDE’s), and
then follow the additional steps below:

(i) Use:
ϑi3 = ϑi3(x

?
1), i = 0, 1, . . . ,m

(ii) For r = 0, 1, . . . ,N − 1, set:

x? = x?r+1 = a+ (r + 1)h, n = 3, k = 3

(iii) For each i = 0, 1, . . . ,m, evaluate:

B{C,3}
i (x?r+1) ≈

n∑
z=0

ϑiz Bk
iz(x

?
r+1)

Output: The approximate solutions for each function:

B{C,3}
0 (x?r+1), B{C,3}

1 (x?r+1), . . . , B{C,3}
m (x?r+1)

III. Third Procedure: Average First and Final Control Point Cubic B-
Spline (FFCPCBS)

After completing the Main Algorithm, proceed with the following steps:

(i) Use:
ϑi3 =

1

2

(
ϑi3(x

?
1) + ϑi3(x

?
N−1)

)
, i = 0, 1, . . . ,m

(ii) For r = 0, 1, . . . ,N − 1, set:

x? = x?r+1 = a+ (r + 1)h, n = 3, k = 3

(iii) For each i = 0, 1, . . . ,m, compute:

B{C,3}
i (x?r+1) ≈

n∑
z=0

ϑiz Bk
iz(x

?
r+1)

Output: The resulting approximations:

B{C,3}
0 (x?r+1), B{C,3}

1 (x?r+1), . . . , B{C,3}
m (x?r+1)



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 19 of 40

6. Numerical Results

In this section, we present some numerical examples to demonstrate the accuracy
and efficiency of the proposed method. We examine a set of test cases to evaluate the
performance of the developed algorithm in solving a system of Volterra integro-differential
equations with classical and fractional orders. The obtained results are compared with
the exact solution; this method’s implementation is supported by a Python program. All
the calculations are executed in least square error.

Example 1. Consider the following classical and fractional order system of Volterra
integro-differential equations (CFVIDEs):

ex Y ′
0(x) + x2 C

aDσ00
x Y0(x) + cos(x)Y0(x) = F0(x) + ω00

∫ x

0
exsY0(s) ds

+ ω01

∫ x

0
(sx2 − 1) CaDβ01

s Y1(s) ds

+ ω02

∫ x

0
s2 sin(x) CaDβ02

s Y2(s) ds,

cos(x)Y ′
1(x) + x3 C

0D0.6
x Y1(x) + ex Y1(x) = F1(x) + ω10

∫ x

0
(x2s+ 1)Y0(s) ds

+ ω11

∫ x

0
s sin(x) C0Dβ11

s Y1(s) ds

+ ω12

∫ x

0
x2(s− 1) C0Dβ12

s Y2(s) ds,

2Y ′
2(x) + ex C

0D0.5
x Y2(x) + sin(x)Y2(x) = F2(x) + ω20

∫ x

0
(x3s2 − 1)Y0(s) ds

+ ω21

∫ x

0
sex C

0Dβ21
s Y1(s) ds

+ ω22

∫ x

0
(s2 + 1) C0Dβ22

s Y2(s) ds.

The given functions F0(x), F1(x), and F2(x) are defined as follows:

F0(x) = ex +
x2.3

Γ(1.3)
+ (x+ 1) cos(x)

− ω00

(
exx3

3
+

exx2

2

)
+

2ω01

Γ(1.2)

(
x4.2

2.2
− x1.2

1.2

)
− 3ω02

Γ(1.5)
· x

3.5

3.5
sin(x),

F1(x) = −2 cos(x)− 2x3.4

Γ(1.4)
+ (1− 2x)ex

− ω10

(
x5

3
+

x4

2
+

x2

2
+ x

)
+

2ω11

Γ(1.3)
· x

2.3

2.3
sin(x)− 3ω12

Γ(1.6)

(
x4.6

2.6
− x3.6

1.6

)
,



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 20 of 40

F2(x) = 6 +
3exx0.5

Γ(1.5)
+ (3x− 1) sin(x)

− ω20

(
x7

4
+

x6

3
− x2

2
− x

)
+

2ω21

Γ(1.1)
· e

xx2.1

2.1
− 3ω22

Γ(1.7)

(
x3.7

3.7
+

x1.7

1.7

)
.

Initial Conditions:

Y0(0) = 1, Y1(0) = 1, Y2(0) = −1

Exact Solutions:

Y0(x) = x+ 1, Y1(x) = 1− 2x, Y2(x) = 3x− 1

Parameters:

σ00 = 0.7, σ10 = 0.6, σ20 = 0.5,

β01 = 0.8, β02 = 0.5, β11 = 0.7,

β12 = 0.4, β21 = 0.9, β22 = 0.3,

ω00 =
sin(0.3)

Γ(13)
, ω01 =

sinh(0.7)

Γ4(16)
, ω02 =

cosh(30)

Γ5(14)
,

ω10 =
cos(89)

Γ(12)
, ω11 =

ln(5)

Γ(12)
, ω12 =

sinh(0.3)

Γ(11)
,

ω20 =
cos(89)

Γ2(9)
, ω21 =

sin(179)

Γ(11)
, ω22 =

sin(30)

Γ(12)
,

N = 10, h = 0.1, x?r = a+ rh, r = 0, 1, . . . ,N − 1.

We aim to approximate the solutions B{C,3}
i (x?r+1), i = 0, 1, 2 as defined in (21). We

execute the programs NFCBI, FCPCBS, and FFCPCBS to compute the unknown control
points ϑi0, ϑi1, ϑi2, ϑi3 for i = 0, 1, 2. These control points are used to construct the
approximate solutions for the given system.

The first table demonstrates values of all control points for B{C,3}
0 (x?r+1), B

{C,3}
1 (x?r+1),

and B{C,3}
2 (x?r+1) according to the three algorithms NFCBI, FCPCBS, and FFCPCBS,

respectively.

Example 2. Consider the following classical and fractional order system of Volterra
integro-differential equations (CFVIDEs):

ex Y ′
0(x) + x4 C

0Dσ00
x Y0(x) + cos(x)Y0(x) = F0(x) + ω00

∫ x

0
(cos(x)− s2)Y0(s) ds

+ ω01

∫ x

0
(3x− s2) C0Dβ01

s Y1(s) ds,

cos(x)Y ′
1(x) + x C

0Dσ10
x Y1(x) + ex Y1(x) = F1(x) + ω10

∫ x

0
(x3 − s4 + 4)Y0(s) ds



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 21 of 40

+ ω11

∫ x

0
(sin(x) + 1) C0Dβ11

s Y1(s) ds.

The given functions F0(x), and F1(x) are defined as follows:

F0(x) = 2xex +
2

Γ(2.3)
x5.3 + (x2 − 1) cos(x)

− ω00

(
x3

3
cos(x)− x cos(x)− x5

5
− x3

3

)
+

2ω01

Γ(2.4)

(
3x3.4

2.4
− x4.4

4.4

)
,

F1(x) = −2x cos(x)− 2

Γ(2.6)
x2.6 + (−x2 + 1)ex

− ω10

(
x6

3
− x7

7
− 4x3

3
− x4 +

x5

5
− 4x

)
− 2ω11

Γ(2.2)

(
x2.2

2.4
sin(x) +

x2.2

2.2

)
.

Initial Conditions:
Y0(0) = −1, Y1(0) = 1

Exact Solutions:
Y0(x) = x2 − 1, Y1(x) = −x2 + 1

Parameters:

σ00 = 0.7, σ10 = 0.4, β01 = 0.6, β11 = 0.8,

ω00 =
sin(0.3)

Γ(13)
, ω01 =

sinh(0.7)

Γ4(16)
, ω10 =

cos(89)

Γ(15)
, ω11 =

ln(2)

Γ(14)
,

N = 10, h = 0.1, x?r = a+ rh, r = 0, 1, . . . ,N − 1.

We aim to approximate the solutions B{C,3}
i (x?r+1), i = 0, 1 as defined in (21). We

execute the programs NFCBI, FCPCBS, and FFCPCBS to compute the unknown con-
trol points ϑi0, ϑi1, ϑi2, ϑi3 for i = 0, 1. These control points are used to construct the
approximate solutions for the given system.

The sixth table demonstrates values of all control points for B{C,3}
0 (x?r+1) and B{C,3}

1 (x?r+1)
according to the three algorithms NFCBI, FCPCBS, and FFCPCBS, respectively.

Example 3. Consider the following classical and fractional order system of Volterra
integro-differential equations (CFVIDEs):

Y ′
0(x) + x C

0D0.4
x Y0(x) + sin(x)Y0(x) = F0(x) + ω00

∫ x

0
xsY0(s) ds

+ ω01

∫ x

0
ex C

0D0.9
s Y1(s) ds,

Y ′
1(x) + x2 C

0D0.5
x Y1(x) + ex Y1(x) = F1(x) + ω10

∫ x

0
(1 + 3x)Y0(s) ds



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 22 of 40

+ ω11

∫ x

0
cos(x) C0D0.7

s Y1(s) ds.

The given functions F0(x) and F1(x) are defined as follows:

F0(x) = 3x2 +
6

Γ(3.6)
x3.6 + (x3 + 1) sin(x)− ω00

(
x6

5
+

x3

2

)
− 2ω01e

x

1.1Γ(1.1)
x1.1

F1(x) = 2 +
2

Γ(1.5)
x2.5 + ex(1 + 2x)− ω10

(
x4

4
+ x+

3x5

4
+ 3x2

)
− 2ω11 cos(x)

1.3Γ(1.3)
x1.3

Initial Conditions:
Y0(0) = 1, Y1(0) = 1

Exact Solutions:
Y0(x) = x3 + 1, Y1(x) = 1 + 2x

Parameters:

σ00 = 0.4, σ10 = 0.5, β01 = 0.9, β11 = 0.7,

ω00 =
sin(0.3)

Γ2(10)
, ω01 =

− ln(0.4× 0.7)

Γ(20)
,

ω10 =
cos(89)

Γ2(10)
, ω11 =

log(0.99999)

Γ(1.5)
,

N = 10, h = 0.1, x?r = a+ rh, r = 0, 1, . . . ,N − 1

We aim to approximate the solutions B{C,3}
i (x?r+1), for i = 0, 1, as defined in equation

(21). We execute the programs NFCBI, FCPCBS, and FFCPCBS to compute the unknown
control points ϑi0, ϑi1, ϑi2, ϑi3 for i = 0, 1. These control points are used to construct the
approximate solutions for the given system.

The 10th table demonstrates values of all control points for B{C,3}
0 (x?r+1) and B{C,3}

1 (x?r+1),
according to the three algorithms: NFCBI, FCPCBS, and FFCPCBS, respectively.

Table 1: The values of control points of B{C,3}
0 (x?

r+1), B
{C,3}
1 (x?

r+1), and B{C,3}
2 (x?

r+1) for three
methods NFCBI, FCPCBS and FFCPCBS, respectively, for example 1.

Method Interval CP of the B{C,3}
0 (x?

r+1) CP of the B{C,3}
1 (x?

r+1) CP of the B{C,3}
2 (x?

r+1)

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.0, 0.1] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9998843842 ϑ13 = −1.000000104 ϑ23 = 2.0000000055

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.1, 0.2] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

Continued on next page



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 23 of 40

Method Interval CP of the B{C,3}
0 (x?

r+1) CP of the B{C,3}
1 (x?

r+1) CP of the B{C,3}
2 (x?

r+1)

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999513774 ϑ13 = −1.000000044 ϑ23 = 2.0000000024

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.2, 0.3] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999738009 ϑ13 = −1.000000024 ϑ23 = 2.0000000014

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.3, 0.4] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999850724 ϑ13 = −1.000000014 ϑ23 = 2.0000000010

NFCBI ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.4, 0.5] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999918765 ϑ13 = −1.000000007 ϑ23 = 2.0000000010

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.5, 0.6] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999964421 ϑ13 = −1.000000003 ϑ23 = 2.0000000011

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.6, 0.7] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999997248 ϑ13 = −1.000000000 ϑ23 = 2.0000000013

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.7, 0.8] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 2.0000022027 ϑ13 = −0.999999997 ϑ23 = 2.0000000017

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.8, 0.9] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 2.0000041419 ϑ13 = −0.999999995 ϑ23 = 2.0000000023

ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.9, 1.0] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996

Continued on next page



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 24 of 40

Method Interval CP of the B{C,3}
0 (x?

r+1) CP of the B{C,3}
1 (x?

r+1) CP of the B{C,3}
2 (x?

r+1)

ϑ03 = 2.000005702 ϑ13 = −0.999999994 ϑ23 = 2.0000000029

FCPCBS ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.0000000000
]0.0, 1.0] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.66667334902 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9998843842 ϑ13 = −1.000000104 ϑ23 = 2.0000000055

FFCPCBS ϑ00 = 1.0000000000 ϑ10 = 1.0000000000 ϑ20 = −1.000000000
]0.0, 1.0] ϑ01 = 1.3333333333 ϑ11 = 0.3333333333 ϑ21 = 0.0000000000

ϑ02 = 1.6666733490 ϑ12 = −0.333333327 ϑ22 = 0.9999999996
ϑ03 = 1.9999450432 ϑ13 = −1.000000049 ϑ23 = 2.0000000042

From (21), we obtain the approximate solution for the classical and fractional order
linear system of Volterra integro-differential equations (SVIDEs-CF) with variable coeffi-
cients, for Example 1 as shown below:

B{C,3}
0 (x?) (NFCBI) =

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

−0.00013566283x?3 + 0.00002004706x?2 + x? + 1 if 0 < x? ≤ 1
10

−0.0000486226x?3 + 0.00002004706x?2 + x? + 1 if 1
10 < x? ≤ 1

5

−0.00002619903x?3 + 0.00002004706x?2 + x? + 1 if 1
5 < x? ≤ 3

10

−0.00001492755x?3 + 0.00002004706x?2 + x? + 1 if 3
10 < x? ≤ 2

5

−0.00000812343x?3 + 0.00002004706x?2 + x? + 1 if 2
5 < x? ≤ 1

2

−0.00000355785x?3 + 0.00002004706x?2 + x? + 1 if 1
2 < x? ≤ 3

5

−0.0000002752x?3 + 0.00002004706x?2 + x? + 1 if 3
5 < x? ≤ 7

10

0.00000220271x?3 + 0.00002004706x?2 + x? + 1 if 7
10 < x? ≤ 4

5

0.00000414190x?3 + 0.00002004706x?2 + x? + 1 if 4
5 < x? ≤ 9

10

0.00000570220x?3 + 0.0000000181x?2 + x? + 1 if 9
10 < x? ≤ 1

B{C,3}
1 (x?) (NFCBI) =

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

−0.0000001223x?3 + 0.00002004706x?2 − 2x? + 1 if 0 < x? ≤ 1
10

−0.0000000442x?3 + 0.00002004706x?2 − 2x? + 1 if 1
10 < x? ≤ 1

5

−0.0000000242x?3 + 0.00002004706x?2 − 2x? + 1 if 1
5 < x? ≤ 3

10

−0.0000000140x?3 + 0.00002004706x?2 − 2x? + 1 if 3
10 < x? ≤ 2

5

−0.0000000077x?3 + 0.00002004706x?2 − 2x? + 1 if 2
5 < x? ≤ 1

2

−0.0000000033x?3 + 0.00002004706x?2 − 2x? + 1 if 1
2 < x? ≤ 3

5

−0.0000000001x?3 + 0.00002004706x?2 − 2x? + 1 if 3
5 < x? ≤ 7

10

0.0000000024x?3 + 0.00002004706x?2 − 2x? + 1 if 7
10 < x? ≤ 4

5

0.0000000042x?3 + 0.00002004706x?2 − 2x? + 1 if 4
5 < x? ≤ 9

10

0.0000000054x?3 + 0.00002004706x?2 − 2x? + 1 if 9
10 < x? ≤ 1



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 25 of 40

B{C,3}
2 (x?) (NFCBI) =

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

0.00000000652x?3 − 0.00000000096x?2 + 3x? − 1 if 0 < x? ≤ 1
10

0.00000000241x?3 − 0.00000000096x?2 + 3x? − 1 if 1
10 < x? ≤ 1

5

0.00000000145x?3 − 0.00000000096x?2 + 3x? − 1 if 1
5 < x? ≤ 3

10

0.00000000108x?3 − 0.00000000096x?2 + 3x? − 1 if 3
10 < x? ≤ 2

5

0.00000000101x?3 − 0.00000000096x?2 + 3x? − 1 if 2
5 < x? ≤ 1

2

0.00000000113x?3 − 0.00000000096x?2 + 3x? − 1 if 1
2 < x? ≤ 3

5

0.00000000139x?3 − 0.00000000096x?2 + 3x? − 1 if 3
5 < x? ≤ 7

10

0.00000000179x?3 − 0.00000000096x?2 + 3x? − 1 if 7
10 < x? ≤ 4

5

0.00000000232x?3 − 0.00000000096x?2 + 3x? − 1 if 4
5 < x? ≤ 9

10

0.00000000298x?3 − 0.00000000096x?2 + 3x? − 1 if 9
10 < x? ≤ 1

FCPCBS =


B{C,3}
0 (x?) = −0.00013566283x?3 + 0.00002004706x?2 + x? + 1 if 0 < x? ≤ 1.0

B{C,3}
1 (x?) = −0.0000001223x?3 + 0.00002004706x?2 − 2x? + 1 if 0 < x? ≤ 1.0

B{C,3}
2 (x?) = 0.00000000652x?3 − 0.00000000096x?2 + 3x? − 1 if 0 < x? ≤ 1.0

FFCPCBS =


B{C,3}
0 (x?) = 1 + x? + 0.00002004706x?2 − 0.00005495675x?3 if 0 < x? ≤ 1.0

B{C,3}
1 (x?) = 1− 2x? + 0.0000000181x?2 − 0.0000000494x?3 if 0 < x? ≤ 1.0

B{C,3}
2 (x?) = −1 + 3x? − 0.00000000096x?2 + 0.00000000427x?3 if 0 < x? ≤ 1.0

Tables 2–4 present a comparison between the approximate and exact solutions of the
functions Y0(x), Y1(x), and Y2(x). The computations are carried out using N = 10, step
size h = 0.1, and sampling points defined by x?r = a + rh, for r = 0, 1, . . . ,N − 1. Three
cubic B-spline methods are considered: NFCBI, FCPCBS, and FFCPCBS, for Example 1
respectively.

Table 2: Compares the exact and approximate based on a least square error of Y0(x) for
Example 1

x?
r Exact B{C,3}

0 (x?
r) (NFCBI) B{C,3}

0 (x?
r) (FCPCBS) B{C,3}

0 (x?
r) (FFCPCBS)

0 1.0 1.0000000000 1.0000000000 1.0000000000
1/10 1.1 1.1000000648 1.1000001689 1.1000001255
1/5 1.2 1.2000002525 1.2000005490 1.2000002019
3/10 1.3 1.3000005556 1.3000009508 1.2999997791
2/5 1.4 1.4000009692 1.4000011846 1.3999984073
1/2 1.5 1.5000014905 1.5000010607 1.4999956363
3/5 1.6 1.6000021183 1.6000003895 1.5999910161
7/10 1.7 1.7000028525 1.6999989813 1.6999840968
4/5 1.8 1.8000036938 1.7999966465 1.7999744282
9/10 1.9 1.9000046433 1.8999931954 1.8999615604
1.0 2.0 2.0000057023 1.9999884384 1.9999450433

L.S.E 8.3881× 10−11 1.9617× 10−10 5.5071× 10−9

Runtime (sec) 2.7607 2.3051 2.2901



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 26 of 40

Table 3: Compares the exact and approximate based on a least square error of Y1(x) for
Example 1

x?
r Exact B{C,3}

1 (x?
r) (NFCBI) B{C,3}

1 (x?
r) (FCPCBS) B{C,3}

1 (x?
r) (FFCPCBS)

0 1.0 1.000000000000 1.000000000000 1.000000000000
1/10 0.8 0.800000000058 0.800000000058 0.800000000113
1/5 0.6 0.600000000222 0.599999999743 0.600000000181
3/10 0.4 0.400000000475 0.399999998322 0.399999999888
2/5 0.2 0.200000000796 0.199999995056 0.199999998564
1/2 0.0 1.16518× 10−9 −1.0782× 10−8 −3.931× 10−9

3/5 -0.2 -0.199999998438 -0.200000019929 -0.200000008090
7/10 -0.4 -0.399999998025 -0.400000033119 -0.400000014319
4/5 -0.6 -0.599999997604 -0.600000051085 -0.600000023023
9/10 -0.8 -0.799999997174 -0.800000074561 -0.800000034606
1.0 -1.0 -0.999999996730 -1.000000104288 -1.000000049477

L.S.E 3.3029× 10−17 2.0681× 10−14 4.4633× 10−15

Runtime (sec) 2.7607 2.3051 2.2901

Table 4: Compares the exact and approximate based on a least square error of Y2(x) for
Example 1

x?
r Exact B{C,3}

2 (x?
r) (NFCBI) B{C,3}

2 (x?
r) (FCPCBS) B{C,3}

2 (x?
r) (FFCPCBS)

0 -1.0 -1.000000000000 -1.000000000000 -1.000000000000
1/10 -0.7 -0.700000000003 -0.700000000003 -0.700000000004
1/5 -0.4 -0.400000000012 -0.399999999986 -0.399999999996
3/10 -0.1 -0.100000000026 -0.099999999911 -0.099999999945
2/5 0.2 0.199999999956 0.200000000265 0.200000000182
1/2 0.5 0.499999999931 0.500000000576 0.500000000415
3/5 0.8 0.799999999903 0.800000001064 0.800000000786
7/10 1.1 1.099999999878 1.100000001777 1.100000001333
4/5 1.4 1.399999999844 1.400000002733 1.400000002077
9/10 1.7 1.699999999889 1.700000003988 1.700000003044
1.0 2.0 1.999999999756 2.000000005566 2.000000004277

L.S.E 1.6268× 10−19 5.8856× 10−17 3.4339× 10−17

Runtime (sec) 2.7607 2.3051 2.2901

Table 5: Least square error for B{C,3}
0 (x?), B{C,3}

1 (x?), and B{C,3}
2 (x?) in Example 1 using different

step sizes.

Method Step size h L.S.E. B{C,3}
0 (x?) L.S.E. B{C,3}

1 (x?) L.S.E. B{C,3}
2 (x?)

NFCBI

1/20 8.3882× 10−13 3.3029× 10−19 1.6271× 10−21

1/50 8.3882× 10−15 3.3028× 10−21 1.6305× 10−23

Continued on next page



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 27 of 40

Method Step size h L.S.E. B{C,3}
0 (x?) L.S.E. B{C,3}

1 (x?) L.S.E. B{C,3}
2 (x?)

1/100 2.0676× 10−15 3.3021× 10−22 1.6875× 10−25

1/1000 2.0677× 10−19 3.2342× 10−27 1.4463× 10−28

FCPCBS

1/20 1.0386× 10−9 8.4106× 10−16 2.4058× 10−18

1/50 6.2009× 10−11 5.0151× 10−17 1.4232× 10−19

1/100 1.8338× 10−11 2.3886× 10−17 6.6535× 10−20

1/1000 4.5237× 10−11 3.0396× 10−17 1.0171× 10−20

FFCPCBS

1/20 2.1963× 10−10 1.4663× 10−16 8.7157× 10−18

1/50 3.3795× 10−11 3.5882× 10−17 4.7505× 10−18

1/100 4.0923× 10−11 4.1016× 10−17 4.6055× 10−18

1/1000 3.8287× 10−11 4.1630× 10−17 4.5875× 10−18

The plots of the approximate solutions obtained by the NFCBI, FCPCBS, and FFCPCBS
methods for Y0(x), Y1(x), and Y2(x) are illustrated in Figures 1, 3, and 5, respectively,
along with the corresponding exact solutions. The least square errors for each iteration
are presented in Figures 2, 4, and 6. These results are based on Example 1 using N = 10
and h = 0.1.

Figure 1: Approximate solution using
NFCBI for Example 1.

Figure 2: Least square error plot using
NFCBI for Example 1.

Figure 3: Approximate solution using
FCPCBS for Example 1.

Figure 4: Least square error plot using
FCPCBS for Example 1.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 28 of 40

Figure 5: Approximate solution using
FFCPCBS for Example 1.

Figure 6: Least square error plot using
FFCPCBS for Example 1.

Table 6: The values of control points of B{C,3}
0 (x?

r+1) and B{C,3}
1 (x?

r+1) for three methods NFCBI,
FCPCBS, and FFCPCBS, respectively, for Example 2

Method Interval CP of B{C,3}
0 (x?

r+1) CP of B{C,3}
1 (x?

r+1)

]0.0, 0.1]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = 1.20000000000011× 10−11

ϑ20 = 1.0000000000000000000
ϑ21 = 1.0000000000000000000
ϑ22 = 0.6666666666608888882
ϑ23 = 1.55986334959834× 10−10

]0.1, 0.2]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = 2.12500000000011× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −5.5622173533720× 10−11

]0.2, 0.3]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = 2.5185185185185765× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −3.35780785892× 10−11

]0.3, 0.4]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666666
ϑ03 = 2.18749999999993× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −2.08461720108132× 10−11

NFCBI ]0.4, 0.5]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = 1.15000000000002× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −3.06465963717528× 10−11

]0.5, 0.6]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = −3.20000000000103× 10−12

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −3.00983517944444× 10−11

Continued on next page



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 29 of 40

Method Interval CP of B{C,3}
0 (x?

r+1) CP of B{C,3}
1 (x?

r+1)

]0.6, 0.7]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = −2.28435374149652× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −3.046413137949× 10−11

]0.7, 0.8]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = −4.572916666667× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −3.08179046518142× 10−11

]0.8, 0.9]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = −8.17622410546180× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −2.96219001164554× 10−11

]0.9, 1.0]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = −1.389125188695603× 10−10

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = −5.600000000000× 10−11

FCPCBS ]0.0, 1.0]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
b03 = 1.20000000000011× 10−11

ϑ20 = 1.00000000000000000000
b21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = 1.559863349598× 10−10

FFCPCBS ]0.0, 1.0]

ϑ00 = −1.0000000000000000000
ϑ01 = −1.0000000000000000000
ϑ02 = −0.6666666666666666667
ϑ03 = −6.345625943477962× 10−11

ϑ20 = 1.00000000000000000000
ϑ21 = 1.00000000000000000000
ϑ22 = 0.66666666666088888882
ϑ23 = 4.999316747991× 10−11

From (21), we obtain the approximate solution for the classical and fractional order
linear system of Volterra integro-differential equations (SVIDEs-CF) with variable coeffi-
cients, for Example 2, as shown below:

B{C,3}
0 (x?) (NFCBI) =

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

1.55986334959834× 10−10x?3 + 0.666666666660882x?2 + 1, 0 < x? ≤ 1
10

−5.56221735337203× 10−11x?3 + 0.666666666660882x?2 + 1, 1
10 < x? ≤ 1

5

−3.35780785892186× 10−11x?3 + 0.666666666660882x?2 + 1, 1
5 < x? ≤ 3

10

−2.08461720108132× 10−11x?3 + 0.666666666660882x?2 + 1, 3
10 < x? ≤ 2

5

−3.06465963717528× 10−11x?3 + 0.666666666660882x?2 + 1, 2
5 < x? ≤ 1

2

−3.00983517944444× 10−11x?3 + 0.666666666660882x?2 + 1, 1
2 < x? ≤ 3

5

−3.04641313794957× 10−11x?3 + 0.666666666660882x?2 + 1, 3
5 < x? ≤ 7

10

−3.08179046518142× 10−11x?3 + 0.666666666660882x?2 + 1, 7
10 < x? ≤ 4

5

−2.96219001164554× 10−11x?3 + 0.666666666660882x?2 + 1, 4
5 < x? ≤ 9

10

−5.60000000000000× 10−11x?3 + 0.666666666660882x?2 + 1, 9
10 < x? ≤ 1



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 30 of 40

B{C,3}
1 (x?) (NFCBI) =

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

1.73339564923936× 10−10x?3 − 1.00000000001735x?2 + 1, 0 < x? ≤ 1
10

−3.82689435696193× 10−11x?3 − 1.00000000001735x?2 + 1, 1
10 < x? ≤ 1

5

−1.6224799281872× 10−12x?3 − 1.00000000001735x?2 + 1, 1
5 < x? ≤ 3

10

−3.49298368007567× 10−10x?3 − 1.00000000001735x?2 + 1, 3
10 < x? ≤ 2

5

−1.32933664076518× 10−10x?3 − 1.00000000001735x?2 + 1, 2
5 < x? ≤ 1

2

−1.27451382780919× 10−10x?3 − 1.00000000001735x?2 + 1, 1
2 < x? ≤ 3

5

−1.31108457424034× 10−10x?3 − 1.00000000001735x?2 + 1, 3
5 < x? ≤ 7

10

−1.34647848426539× 10−10x?3 − 1.00000000001735x?2 + 1, 7
10 < x? ≤ 4

5

−1.22686305559228× 10−10x?3 − 1.00000000001735x?2 + 1, 4
5 < x? ≤ 9

10

−3.86468634872017× 10−10x?3 − 1.00000000001735x?2 + 1, 9
10 < x? ≤ 1

B{C,3}
0 (x?) (FCPCBS) = 1.199929045014870× 10−11x?3 + 1.000000x?2 − 1, 0 < x? ≤ 1

B{C,3}
1 (x?) (FCPCBS) = 1.73339564923936× 10−10x?3 − 1.00000000001735x?2 + 1, 0 < x? ≤ 1

B{C,3}
0 (x?) (FFCPCBS) = −6.34567953738951× 10−11x?3 + 1.000000x?2 − 1, 0 < x? ≤ 1

B{C,3}
1 (x?) (FFCPCBS) = 6.7346128673762× 10−11x?3 − 1.00000000001735x?2 + 1, 0 < x? ≤ 1

Table 7: Compares the exact and approximate based on a least square error of Y0(x) for
Example 2.

x?
r Exact B{C,3}

0 (x?
r) (NFCBI) B{C,3}

0 (x?
r) (FCPCBS) B{C,3}

0 (x?
r) (FFCPCBS)

0 -1.00 -1.00 -1.00 -1.00
1/10 -0.99 -0.99 -0.99 -0.99
1/5 -0.96 -0.96 -0.96 -0.960000000001
3/10 -0.91 -0.909999999999 -0.910000000000 -0.910000000002
2/5 -0.84 -0.839999999998 -0.839999999999 -0.840000000004
1/2 -0.75 -0.749999999995 -0.749999999998 -0.750000000008
3/5 -0.64 -0.639999999992 -0.639999999997 -0.640000000014
7/10 -0.51 -0.509999999986 -0.509999999996 -0.510000000022
4/5 -0.36 -0.359999999980 -0.359999999994 -0.360000000032
9/10 -0.19 -0.189999999972 -0.189999999991 -0.190000000046
1.0 0.00 3.61087481×10−11 1.2000000×10−11 -6.3000000×10−11

L.S.E 8.3881× 10−11 2.7712× 10−21 2.8489× 10−22

Run Time (s) 2.7607 0.9465 1.1776



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 31 of 40

Table 8: Compares the exact and approximate based on a least square error of Y1(x) for
Example 2.

x?
r Exact B{C,3}

1 (x?
r) (NFCBI) B{C,3}

1 (x?
r) (FCPCBS) B{C,3}

1 (x?
r) (FFCPCBS)

0 1.00 1.00 1.00 1.00
1/10 0.99 0.99 0.99 0.99
1/5 0.96 0.959999999999 0.960000000001 0.96
3/10 0.91 0.909999999998 0.910000000003 0.91
2/5 0.84 0.839999999997 0.840000000008 0.840000000002
1/2 0.75 0.749999999994 0.750000000017 0.750000000004
3/5 0.64 0.639999999991 0.640000000031 0.640000000008
7/10 0.51 0.509999999987 0.510000000051 0.510000000015
4/5 0.36 0.359999999982 0.360000000078 0.360000000023
9/10 0.19 0.189999999977 0.190000000112 0.190000000035
1.0 0.00 -2.8274724×10−11 1.560000×10−10 5.000000×10−11

L.S.E 2.7712× 10−21 1.9151× 10−21 4.6922× 10−21

Run Time (s) 0.9465 0.9465 1.1776

Table 9: Least square error for Y0(x), and Y1(x), in Example 2 using different step sizes.

Method Step Size h LSE of B{C,3}
0 (x?) LSE of B{C,3}

1 (x?)

NFCBI

1/10 2.7712× 10−21 1.9151× 10−21

1/40 2.1379× 10−23 2.1333× 10−23

1/80 5.8327× 10−25 1.1559× 10−24

1/100 1.4706× 10−25 1.0112× 10−24

FCPCBS

1/10 2.8489× 10−22 4.6922× 10−21

1/40 2.8487× 10−24 4.8138× 10−22

1/80 2.8303× 10−28 4.8167× 10−26

1/100 2.8490× 10−30 2.2284× 10−28

FFCPCBS

1/10 7.9664× 10−21 4.5745× 10−21

1/40 3.7361× 10−21 1.7798× 10−21

1/80 5.6535× 10−23 3.1644× 10−23

1/100 3.5825× 10−24 7.9086× 10−24

The plots of the approximate solutions obtained by NFCBI, FCPCBS, and FFCPCBS for
Y0(x) and Y1(x) are shown in Figures 7, 9, and 11, respectively, along with the exact solution
plots. The corresponding least square errors for each iteration are illustrated in Figures 8, 10, and
12, respectively, for Example 2.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 32 of 40

Figure 7: Approximation solution of Y0(x)
and Y1(x) method NFCBI for example 2.

Figure 8: Least square error plot of Y0(x)
and Y1(x) method NFCBI for example 2.

Figure 9: Approximation solution of Y0(x)
and Y1(x) method FCPCBS for example 2.

Figure 10: Least square error plot of Y0(x)
and Y1(x) method FCPCBSfor example

6.2.

Figure 11: Approximation solution of Y0(x)
and Y1(x) method FFCPCBS for

example 2.

Figure 12: Least square error plot of Y0(x)
and Y1(x) method FFCPCBS for example

2.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 33 of 40

Table 10: The values of control points of B{C,3}
0 (x?

r+1) and B{C,3}
1 (x?

r+1) for three methods NFCBI,
FCPCBS, and FFCPCBS, respectively, for Example 3.

Method Interval CP of B{C,3}
0 (x?

r+1) CP of B{C,3}
1 (x?

r+1)

]0.0, 0.1]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999991827737578

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999988454

]0.1, 0.2]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999996593810891

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999992006

]0.2, 0.3]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999998175375771

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999995115

]0.3, 0.4]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999998962879832

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999994227

NFCBI ]0.4, 0.5]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999999433926590

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999989342

]0.5, 0.6]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999999747462898

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999985345

]0.6, 0.7]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999999971441733

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999981348

]0.7, 0.8]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 2.000000000139726009

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999974687

Continued on next page



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 34 of 40

Method Interval CP of B{C,3}
0 (x?

r+1) CP of B{C,3}
1 (x?

r+1)

]0.8, 0.9]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 2.000000000271051181

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999969358

]0.9, 1.0]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 2.000000000376601861

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999964917

FCPCBS ]0.0, 1.0]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999991827737578

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999988454

FFCPCBS ]0.0, 1.0]

ϑ00 = 1.000000000000000000
ϑ01 = 1.000000000000000000
ϑ02 = 1.000000000478453520
ϑ03 = 1.999999996102169719

ϑ20 = 1.000000000000000000
ϑ21 = 1.666666666666666666
ϑ22 = 2.333333333333333335
ϑ23 = 2.999999999999976907

From (21), we obtain the approximate solution for the classical and fractional order linear sys-
tem of Volterra integro-differential equations (SVIDEs-CF) with variable coefficients, for Example
3, as shown below:

B{C,3}
0 (x?) (NFCBI) =

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

0.999999990392377x?3 + 1.43536027508162× 10−9 x?2 + 1, 0 < x? ≤ 1
10

0.999999995158451x?3 + 1.43536027508162× 10−9 x?2 + 1, 1
10 < x? ≤ 1

5

0.999999996740016x?3 + 1.43536027508162× 10−9 x?2 + 1, 1
5 < x? ≤ 3

10

0.999999997527520x?3 + 1.43536027508162× 10−9 x?2 + 1, 3
10 < x? ≤ 2

5

0.999999997998566x?3 + 1.43536027508162× 10−9 x?2 + 1, 2
5 < x? ≤ 1

2

0.999999998312103x?3 + 1.43536027508162× 10−9 x?2 + 1, 1
2 < x? ≤ 3

5

0.999999998536081x?3 + 1.43536027508162× 10−9 x?2 + 1, 3
5 < x? ≤ 7

10

0.999999998704366x?3 + 1.43536027508162× 10−9 x?2 + 1, 7
10 < x? ≤ 4

5

0.999999998835690x?3 + 1.43536027508162× 10−9 x?2 + 1, 4
5 < x? ≤ 9

10

0.999999998941242x?3 + 1.43536027508162× 10−9 x?2 + 1, 9
10 < x? ≤ 1



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 35 of 40

B{C,3}
1 (x?) (NFCBI) =

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

−1.15463194561016× 10−14 x?3 + 2.00000000000000x? + 1, 0 < x? ≤ 1
10

−7.99360577730113× 10−15 x?3 + 2.00000000000000x? + 1, 1
10 < x? ≤ 1

5

−5.32907051820075× 10−15 x?3 + 2.00000000000000x? + 1, 1
5 < x? ≤ 3

10

−5.32907051820075× 10−15 x?3 + 2.00000000000000x? + 1, 3
10 < x? ≤ 2

5

−1.06581410364015× 10−14 x?3 + 2.00000000000000x? + 1, 2
5 < x? ≤ 1

2

−1.42108547152020× 10−14 x?3 + 2.00000000000000x? + 1, 1
2 < x? ≤ 3

5

−1.86517468137026× 10−14 x?3 + 2.00000000000000x? + 1, 3
5 < x? ≤ 7

10

−2.48689957516035× 10−14 x?3 + 2.00000000000000x? + 1, 7
10 < x? ≤ 4

5

−3.01980662698043× 10−14 x?3 + 2.00000000000000x? + 1, 4
5 < x? ≤ 9

10

−3.55271367880050× 10−14 x?3 + 2.00000000000000x? + 1, 9
10 < x? ≤ 1

B{C,3}
0 (x?) (FCPCBS) = 0.999999990392377x?3 + 1.43536027508162× 10−9 x?2 + 1, 0 < x? ≤ 1.0

B{C,3}
1 (x?) (FCPCBS) = −1.15463194561016× 10−14 x?3 + 2.00000000000000x? + 1, 0 < x? ≤ 1.0

B{C,3}
0 (x?) (FFCPCBS) = 0.999999994666809x?3 + 1.43536027508162× 10−9 x?2 + 1, 0 < x? ≤ 1.0

B{C,3}
1 (x?) (FFCPCBS) = −2.30926389122033× 10−14 x?3 + 2.00000000000000x? + 1, 0 < x? ≤ 1.0

Tables 11 and 12 demonstrate a comparison between the approximate and exact solutions of
Y0(x) and Y1(x), respectively. The computations are carried out by setting N = 10, h = 0.1, and
x?
r = a + rh for r = 0, 1, . . . ,N − 1, using three cubic B-spline methods: NFCBI, FCPCBS, and

FFCPCBS.

Table 11: Comparison of the exact and approximate solutions based on least square error
for Y0(x) for Example 3.

x?
r Exact B{C,3}

0 (x?
r) (NFCBI) B{C,3}

0 (x?
r) (FCPCBS) B{C,3}

0 (x?
r) (FFCPCBS)

0 1.000 1.000 1.000 1.000
1/10 1.001 1.00100000000475 1.001000000005 1.001000000009
1/5 1.008 1.00800000001868 1.007999999981 1.008000000015
3/10 1.027 1.02700000004116 1.026999999870 1.026999999985
2/5 1.064 1.06400000007142 1.063999999615 1.063999999888
1/2 1.125 1.12500000010866 1.124999999158 1.124999999692
3/5 1.216 1.21600000015215 1.215999998441 1.215999999365
7/10 1.343 1.34300000020120 1.342999997408 1.342999998874
4/5 1.512 1.51200000025526 1.511999996000 1.511999998188
9/10 1.729 1.72900000031386 1.728999994159 1.728999997275
1.0 2.000 2.00000000037660 1.999999991828 1.999999996102

L.S.E 3.881× 10−19 1.269× 10−16 2.768× 10−17

Runtime (s) 1.2066 1.2065 1.2085



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 36 of 40

Table 12: Comparison of the exact and approximate solutions based on least square error for
Y1(x) for Example 3.

x?
r Exact B{C,3}

1 (x?
r) (NFCBI) B{C,3}

1 (x?
r) (FCPCBS) B{C,3}

1 (x?
r) (FFCPCBS)

0 1.0 1.0 1.0 1.0
1/10 1.2 1.2 1.2 1.2
1/5 1.4 1.4 1.4 1.4
3/10 1.6 1.6 1.6 1.6
2/5 1.8 1.8 1.8 1.80000000000009
1/2 2.0 2.0 2.0000000000008 2.0000000000007
3/5 2.2 2.2 2.2000000000007 2.2000000000017
7/10 2.4 2.4 2.4000000000099 2.4000000000058
4/5 2.6 2.6 2.6000000000023 2.6000000000014
9/10 2.8 2.79999999999999 2.7999999999998 2.7000000000238
1.0 3.0 3.00000000000001 3.00000000000013 3.00000000000013

L.S.E 8.239× 10−29 2.598× 10−28 1.060× 10−27

Runtime (s) 1.2066 1.2065 1.2085

Table 13: Least square error for Y0(x), and Y1(x), in Example 3 using different step sizes.

Method Step Size h LSE of B{C,3}
0 (x?) LSE of B{C,3}

1 (x?)

NFCBI
1/10 3.88101892329981× 10−19 8.23866607890194× 10−29

1/50 2.43602430852715× 10−22 7.13328499648335× 10−29

1/1000 1.38050658413677× 10−29 2.23866607890190× 10−29

FCPCBS
1/10 1.2693341933429862× 10−16 2.5983106065717076× 10−28

1/50 7.966390404945247× 10−20 1.2232274411583314× 10−28

1/1000 3.944304526105059× 10−31 3.4512664603419266× 10−31

FFCPCBS
1/10 2.768234509326171× 10−17 1.0601304490038872× 10−27

1/50 6.488786754097752× 10−19 1.0597574191466658× 10−27

1/1000 1.7749370367472766× 10−29 1.0601304490038872× 10−28

The plot of these approximate solutions obtained by NFCBI, FCPCBS and FFCPCBS for Y0(x)
, and Y1(x) are demonstrate in figures 13,15 and 17 with the exact solution plots, respectively.
The least square error of each iteration are shows in figures 14, 16 and 18 for example 3.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 37 of 40

Figure 13: Approximation solution of Y0(x)
and Y1(x) method NFCBI for example 3.

Figure 14: Least square error plot of Y0(x)
and Y1(x) method NFCBI for example 5.3.

Figure 15: Approximation solution of Y0(x)
and Y1(x) method FCPCBS for example 3.

Figure 16: Least square error plot of Y0(x)
and Y1(x) method FCPCBSfor example 3.

Figure 17: Approximation solution of
Y0(x)and Y1(x) method FFCPCBS for

example 3.

Figure 18: Least square error plot of Y0(x)
and Y1(x) method FFCPCBS for example

3.

7. Conclusion

In this work, a high-accuracy cubic B-spline collocation method is presented for solving a
system of Volterra integro-differential equations with classical and fractional Caputo derivatives
(SVIDE’s-CF). Unlike previous studies that primarily focus on single equations [29, 32–34], the



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 38 of 40

proposed approach efficiently handles systems of equations. The method was implemented in
Python and validated through several benchmark examples (Tables 5, 9, and 13). The resulting
least-squares errors demonstrate improvements by several orders of magnitude, highlighting the
method’s superior numerical accuracy and stable convergence properties across a variety of test
problems. To further develop the method, future research can explore its generalization to non-
linear systems, time-dependent parameters, and non-local operators. Additionally, incorporating
adaptive mesh generation and refinement algorithms can render it more versatile and efficient. In-
deed, this method is a powerful computational technique with applications across many disciplines
of applied mathematics, physics, and engineering.

Acknowledgements

This work is part of the PhD research conducted by the student Diar Khalid Abdullah in the
Department of Mathematics, College of Science, University of Sulaimani, Iraq.

Conflict of Interest

The authors declare no conflict of interest.

Availability of data material

All the data of the study are included.

Authors’ Contributions

D.Kh.A., Sh.Sh.A. and K. H.F. wrote the main manuscript. All authors reviewed the manuscript.

References

[1] A. Domoshnitsky, A. Sitkin, and L. Zuckerman. Mathematical modeling of covid-19 trans-
mission in the form of system of integro-differential equations. Mathematics, 10(23):4500,
2022.

[2] A. Turab, A. Montoyo, and J. A. Nescolarde-Selva. Computational and analytical analysis of
integral-differential equations for modeling avoidance learning behavior. Journal of Applied
Mathematics and Computing, 70(5):4423–4439, 2024.

[3] I. Podlubny. Fractional Differential Equations, volume 198 of Mathematics in Science and
Engineering. Academic Press, USA, 1999.

[4] H. Mohammadi, S. Kumar, Sh. Rezapour, and S. Etemad. A theoretical study of the caputo-
fabrizio fractional modeling for hearing loss due to mumps virus with optimal control. Chaos,
Solitons & Fractals, 144:110668, 2021.

[5] S. Alshammari, M. Alshammari, M. Alabedalhadi, M. M. Al-Sawalha, and M. Al-Smadi.
Numerical investigation of a fractional model of a tumor-immune surveillance via caputo
operator. Alexandria Engineering Journal, 86:525–536, 2024.

[6] I. Ali, M. Yaseen, and I. Akram. Utilizing cubic b-spline collocation technique for solving
linear and nonlinear fractional integro-differential equations of volterra and fredholm types.
Fractal and Fractional, 8(5):268, 2024.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 39 of 40

[7] K. Hama Faraj, Sh. Shawki, and D. Khalid. Numerical treatment solution of volterra integro-
fractional differential equation by using linear spline function. Journal of Zankoy Sulaimani
- Part A, 22(2), 2020.

[8] J. Chang, Q. Yang, and C. Liu. B-spline method for solving boundary value problems of
linear ordinary differential equations. In Advances in Computer Science and Information
Engineering, pages 326–333. Springer, 2010.

[9] M. Yi and J. Huang. Cas wavelet method for solving the fractional integro-differential
equation with a weakly singular kernel. International Journal of Computer Mathematics,
92(8):1715–1728, 2015.

[10] K. R. Raslan and T. S. El-Danaf. Solitary waves solutions of the mrlw equation using quintic
b-splines. Journal of King Saud University - Science, 22(3):161–166, 2010.

[11] F. Haq, S. Islam, and I. Tirmizi. A numerical technique for solution of the mrlw equation
using quartic b-splines. Applied Mathematical Modelling, 34(12):4151–4160, 2010.

[12] K. Hama Faraj, Sh. Shawki, and D. Khalid. Approximate solution of volterra integro-fractional
differential equations using quadratic spline function. Bulletin of Karaganda University. Math-
ematics Series, 101(1):50–64, 2021.

[13] M. Ebrahimi and S. Rouhani Rankoohi. Spline collocation for system of fredholm and volterra
integro-differential equations. Knowledge and Information Systems, 43(1):219–247, 2015.

[14] M. Abbas, S. Aslam, F. A. Abdullah, M. B. Riaz, and K. A. Gepreel. An efficient spline
technique for solving time-fractional integro-differential equations. Heliyon, 9(9), 2023.

[15] B. Ghosh and J. Mohapatra. An iterative difference scheme for solving arbitrary order non-
linear volterra integro-differential population growth model. Journal of Analysis, 2023.

[16] I. Masti and K. Sayevand. On collocation-galerkin method and fractional b-spline functions
for a class of stochastic fractional integro-differential equations. Mathematics and Computers
in Simulation, 216:263–287, 2024.

[17] Q. Khan, H. Khan, P. Kumam, F. Tchier, and G. Singh. Ladm procedure to find the ana-
lytical solutions of the nonlinear fractional dynamics of partial integro-differential equations.
Demonstratio Mathematica, 57(1), 2024.

[18] M. Khader and N. Sweilam. On the approximate solutions for system of fractional integro-
differential equations using chebyshev pseudo-spectral method. Applied Mathematical Mod-
elling, 37(24):9819–9828, 2013.

[19] L. Huang, X. Li, Y. Zhao, and X. Duan. Approximate solution of fractional integro-
differential equations by taylor expansion method. Computers & Mathematics with Appli-
cations, 62(3):1127–1134, 2011.

[20] H. Yu and W. Zhang. Solving fractional integro-differential equations with singular kernels
using the fractional taylor series method. Journal of Computational and Applied Mathematics,
421, 2023.

[21] A. Singh, R. Gupta, and T. Patel. An adaptive mesh method for fractional integro-differential
equations. Numerical Mathematics Methods and Applications, 17(1):55–78, 2024.

[22] L. Li and Y. Liu. Fractional integro-differential equations: Numerical solutions using deep
learning techniques. Journal of Computational Physics, 2024.

[23] Sh. Shawki. On System of Linear Volterra Integro-Fractional Differential Equations. Sulaimani
University, Sulaimani, 2006.

[24] K. Mariya. Properties and applications of the Caputo fractional operator. Karlsruhe, Bulgaria,
2005.

[25] B. Rajani and D. Debasish. A mixed quadrature rule by blending clenshaw-curtis and gauss-
legendre quadrature rules for approximation of real definite integrals in adaptive environment.
International Association of Engineers, 1, 2011.



D. Kh. Abdullah, Sh. Sh. Ahmed, K. H. F. Jwamer / Eur. J. Pure Appl. Math, 18 (4) (2025), 6763 40 of 40

[26] M. Sun and Q. Wu. On the chebyshev spectral collocation method for the solution of highly
oscillatory volterra integral equations of the second kind. Applied Mathematics and Nonlinear
Sciences, 9(1), 2024.

[27] N. Ronald and L. Tom. Knot insertion and deletion algorithms for B-spline curves and
surfaces. Library of Congress Cataloging in Publication, 1993.

[28] Z. Mahmoodi, J. Rashidinia, and E. Babolian. B-spline collocation method for linear and non-
linear fredholm and volterra integro-differential equations. Applied Analysis, 92(9):1787–1802,
2013.

[29] M. Gholamian and J. Saberi-Nadjafi. Cubic b-splines collocation method for a class of partial
integro-differential equation. Alexandria Engineering Journal, 57(3):2157–2165, 2018.

[30] B. Sambhunath and C. Brain. Bézier and Splines in Image Processing and Machine Vision.
Springer, London, 2008.

[31] D. Khalid, Sh. Shawki, and K. Hama Faraj. Numerical approximation of solving volterra
integro-fractional differential equations using b-spline functions. Journal of Southwest Jiaotong
University, 59(4), 2024.

[32] S. Khan and M. Misro. Numerical applications of cubic b-spline collocation methods: A
systematic review. page 030019, 2024.

[33] M. Farshid and A. Sahar. Cubic b-spline approximation for linear stochastic integro-
differential equation of fractional order. Journal of Computational and Applied Mathematics,
366, 2020.

[34] F. Mirzaee and S. Alipour. Cubic b-spline approximation for linear stochastic integro-
differential equation of fractional order. Journal of Computational and Applied Mathematics,
336:112440, 2020.