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

A Method Employing Legendre Wavelets and a Finite Iterative
Approach for Efficiently Solving Systems of Linear Fredholm

Integral Equations
Mohamed A. Ramadan1, Mohamed Adel2,∗, Heba M. Arafa3

1 Department of Mathematics and Computer Science, Faculty of Science, Menoufia University, Menoufia 32511,
Shebein El Kom, Egypt
2 Department of Mathematics, Faculty of Science, Islamic University of Madinah, 42210, Medina, KSA
3 Department of Mathematics, Faculty of Education, Ain Shams University, Roxy 11341, Cairo, Egypt

Abstract. Integral equations are essential in numerous domains of practical mathematics. This article introduces
a straightforward, precise, and efficient iterative technique for resolving one-dimensional Fredholm integral equations
of the second class. The suggested numerical method relies on Legendre wavelet functions. Utilizing these wavelets,
the integral equation system is converted into a duo of interconnected systems of algebraic matrix equations. A finite
iterative approach is employed to resolve these systems and ascertain the coefficients that formulate the approximate
numerical solutions of the unknown functions. A variety of examples are included to evaluate the proposed numerical
approach.
2020 Mathematics Subject Classifications: 45B05, 65R20, 65T60, 65Y20, 42C40
Key Words and Phrases: Fredholm integral equations of the second kind, one-dimensional integral systems, Leg-
endre wavelets, iterative methods, numerical approximation, matrix equations, wavelet-based methods

1. Introduction

Integral equations are essential in the mathematical modeling of numerous applied science and engineer-
ing challenges, encompassing areas such as quantum physics, thermal analysis, and signal processing ([1–3]).
Due to the frequent inaccessibility of precise solutions for several integral equations, researchers have devel-
oped a diverse array of analytical and computational methods to approximate solutions for various types
of integral equations ([4–11]). These techniques often include orthogonal polynomials, including Legendre,
Bernstein, Jacobi, Chebyshev, and triangular functions, and Laguerre polynomials, which are extensively
applied in resolving integrals that involve special functions.

Integral equation systems (IES), especially of the Fredholm variety, typically cannot be resolved in closed
form, requiring effective numerical methods. Researchers have suggested a diverse range of such method-
ologies. Babolian et al. ([12, 13]) presented direct and decomposition methods; Jafarian and colleagues
([14, 15]) utilized neural networks and Bernstein collocation; and Mahmoodi [16] implemented collocation
and spectral methods for both Fredholm and Volterra systems. Alipour et al. [17] expanded these concepts
to coupled systems, whereas Huang et al. [18] employed Taylor expansion. Maleknejad et al. [19] utilized
block pulse functions (BPFs). Ramadan et al. ([20, 21]) proposed extended and generalized finite iterative
solution methodologies. Maleknejad et al. [22] achieved novel results utilizing Taylor series for first-kind sys-
tems. In [23], Ramadan et al. expanded the triangle function in conjunction with a finite iterative technique.
Maleknejad et al. [24] employed second-kind Taylor series systems. H. Almasieh and Roodaki [25] employed

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

Email addresses: mohamed.Abdellatif@science.menofia.edu.eg (M. A. Ramadan),
adel@sci.cu.edu.eg (M. Adel), hebaallahmohammed@edu.asu.edu.eg (H. M. Arafa)

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



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 2 of 18

triangular functions to solve a system of Fredholm integral equations, whereas Golbabai and Keramati [26]
proposed an effective computational methodology. Alternative methodologies encompass triangle function
techniques and B-spline wavelet frameworks ([27, 28]). Najafi et al. [29] introduced a linear Legendre multi-
wavelets method for addressing systems of Fredholm integral equations, while Z. Elahiet et al. [30] applied
the Laguerre method for solving linear systems of Fredholm integral equations. Additionally, Elahi et al.
([31, 32]) utilized Laguerre and Bessel polynomials for integro-differential and differential-difference equa-
tions. Recently, Ramadan et al. [31] introduced a method utilizing triangular functions to solve systems of
linear Fredholm integral equations by an efficient finite iterative algorithm, while Arafa and Ramadan [33]
and Legendre wavelet-based methods were presented for addressing coupled Fredholm systems.
The structure of this paper is as follows:Section 1 provides the introduction. In Section 2, we present key
definitions and preliminaries, including an overview of Legendre wavelets, their properties, function approx-
imation, and an iterative algorithm for solving coupled system matrix equations. Section 3 outlines our
proposed method for addressing linear systems of Fredholm integral equations of the second kind. We pro-
vide the convergence analysis and error estimation in Section 4. Section 5 presents three numerical examples
that showcase the precision and effectiveness of the proposed method. Finally, Section 5 concludes the paper
by summarizing the main findings.

2. Definitions and Preliminaries

2.1. Legendre wavelet and its properties

Considering a single function ”mother wavelet” ψ(t), from which wavelets represent a family of functions
by dilating and transforming this single function. This family of continuous wavelets [17] has the following
form:

ψa,b(t) = |a|−
1
2

(
t− b

a

)
, a, b ∈ R, a 6= 0 (2.1)

The Legendre wavelets on the interval [ 0,1 ) defined by

ψn,m(t) =

{√
m+ 1

22
k
2Lm

(
2kt− n̂

)
, n̂−1

2k
≤ t < n̂+1

2k
,

0, otherwise,
(2.2)

for which k is positive integer, n = 1, 2, . . . , 2k−1 and n̂ = 2n− 1, the order of the Legendre Polynomial
is denoted by m = 0, 1, 2, . . . ,M − 1 and the normalized time is denoted by t. The Legendre Polynomials
Lm which are obtained in the above definition is proposed as follows:

L0(t) = 1,

L1(t) = t,

Lm+1(t) =
2m+ 1

m+ 1
tLm(t)− m

m+ 1
Lm−1(t) , m = 1, 2, 3, . . . (2.3)

which are orthogonal over [-1,1] with weighting function.

2.2. Function Approximation

A function f(t) which is defined on [0, 1) can be extended as Legendre Wavelet infinite series of the
following type

f(t) =
∞∑
n=1

∞∑
m=0

cn,mψn,m, (2.4)



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 3 of 18

where cn,m = 〈f, ψn,m〉.
After being trimmed, Eq. (2.4) can be rewritten as follows.

f(t) ≈
2k−1∑
n=1

cn,m

M∑
m=1

ψn,m = CTψ(t), (2.5)

where C =
[
cl,0 cl,1 . . . cl,M−1 . . . c2k−1,0 c2k−1,1 · · · c2k−1,M−1

]T
,

and
ψ(t) =

[
ψ1,0 ψ1,1 . . . ψ1,M−1 . . . ψ2k−1,0 ψ2k−1,1 . . . ψ2k−1,M−1

]T .

2.3. Iterative algorithm for solving the coupled system matrix equations

Finite iterative methods are a class of numerical algorithms designed to solve matrix equations and
coupled matrix equations by reaching the exact solution in a finite number of steps, assuming exact arith-
metic. Unlike traditional iterative methods that converge asymptotically, these methods terminate after a
predetermined number of iterations, making them highly efficient and predictable in computational cost. For
example, see Ramadan et al. [20] Extended and generalized finite iterative solution techniques previously
applied to specific Sylvester-type equations. Moreover, Ramdan et al. [21] developed an iterative method
tailored for bisymmetric structured solutions, focused on least-norm generalized solutions. Ramadan and
Ali [10] used orthogonal triangular functions to convert the linear system of Fredholm integral equations into
a coupled matrix equation, then applies a finite iterative algorithm. In addition, Ramadan et al. [23] ex-
tended the triangular function coupled with finite iterative method to 2D fuzzy Fredholm integral equations,
transforming them into a coupled algebraic system and solving with finite iteration. Recently, Ramadan et
al. [33] approximated the solution of system of linear Fredholm integral equations via a set of orthogonal
triangular functions, which converts the continuous integral equations into four coupled algebraic matrix
equations. These matrix equations are then tackled with an efficient finite-step iterative algorithm. Based
on the finite-step iterative algorithms developed in the literature, particularly those introduced by M.A. Ra-
madan et al. [10,20,2123,33] and other related works, we propose in this subsection a new iterative algorithm
for solving a coupled system of linear matrix equations of the forms:

A1C1 +B1C2 = F and A2C1 +B2C2 = G.

Algorithm 2.1
1- Input A1, B1, A2, B2, F,G
2- Choose arbitrary vectors C11 ∈ Cn×1 and C21 ∈ Cn×1

3- Set

R1 = diag (F − f (C11 , C21) , (G− g (C11 , C21)) ,

S1 = AT
1 ( F− f (C11 , C21)) +AT

2 (G− g (C11 , C21)) ,

T1 = BT
1 ( F− f (C11 , C21)) +BT

2 (G− g (C11 , C21)) ,

where

f (C11 , C11) = A1C11 +B1C21

g (C11 , C21) = A2C11 +B2C21

4- If Rk = 0 then stop and C1k , Y C2k is the solution else let k = k + 1 go to step 5 .
5- Compute



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 4 of 18

C1k+1
= C1k +

‖Rk‖2

‖Sk‖2 + ‖Tk‖2
Sk, C2k+1

= C2k +
‖Rk‖2

‖Sk‖2 + ‖Tk‖2
Tk,

Rk+1 = diag
(
F− f

(
C1k+1

, C2k+1

)
,
(
G− g

(
C1k+1

, C2k+1

))
= Rk −

‖Rk‖2

‖Sk‖2 + ‖Tk‖2
diag (f (Sk, Tk) , g (Sk, Tk)) ,

Sk+1 = AT
1

(
F − f

(
Y1k+1

, Y2k+1

))
+AT

2

(
G− g

(
C1k+1

, C2k+1

))
+

‖Rk+1‖2

‖Rk‖2
Sk,

Tk+1 = V
(
F− f

(
Y1k+1

, Y2k+1

))
+BT

2

(
G− g

(
C1k+1

, C2k+1

))
+

‖Rk+1‖2

‖Rk‖2
Tk.

3. Legendre wavelet method for linear Fredholm integral equations system of second
kind

The proposed approach begins by using the Legendre wavelet method to convert the coupled linear
system if Fredholm integral equations into a system of coupled matrix algebraic equations. Next, Algorithm
2.1 is employed to solve the resulting system of Sylvester matrix equations, yielding the solution function
for the original problem. First, consider the following equations:

u(x) = f(x) + λ1

(∫ 1

0
(k1(x, t)u(t)) dt+

∫ 1

0
(k2(x, t)v(t))

)
dt

v(x) = g(x) + λ2

(∫ 1

0
(k3(x, t)u(t)) dt+

∫ 1

0
(k4(x, t)v(t))

)
dt (3.1)

where k1(x, t), k2(x, t) ∈ L2([0, 1]× [0, 1]) and h(x), g(x) ∈ L2([0, 1])
The unknown functions u(x), v(x) can be expanded as

u(x) ≈ CT
1 ψ(t), v(x) ≈ CT

2 ψ(t) (3.2)

where C1 and C2 are the unknown 2k−1M vectors and ψ(t) is given by Eq. (2.5) and (2.6). Likewise,
k1(x, t), k2(x, t), k3(x, t), k4(x, t), f(x), g(x) are also expanded into the LWM as:

k1(x, t) ≈ ψT(x)K1ψ(t), k2(x, t) ≈ ψT(x)K2ψ(t)

k3(x, t) ≈ ψT(x)K3ψ(t), k4(x, t) ≈ ψT(x)K4ψ(t)

f(x) ≈ F Tψ(x), g(x) ≈ GTψ(x) (3.3)

After substituting the approximate equations (3.2) - (3.3) into (3.1) we get,

ψT (x)C1 = ψT (x)F + λ1ψ
T (x)K1

∫ 1

0
ψ(t)ψT (t)C1dt

+λ1∅T (x)K2

∫ 1

0
ψ(t)ψT (t)C2dt

ψT (x)C2 = ψT (x)G+ λ2ψ
T (x)K3

∫ 1

0
ψ(t)ψT (t)C1dt (3.4)

+λ2ψ
T (x)K4

∫ 1

0
ψ(t)ψT (t)C2dt



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 5 of 18

where, ∫ 1

0
ψ(t)ψT(t)dt = I (3.5)

Making use of Eq. (3.5), we get,

ψT(x)C1 = F Tψ(x) + +λ1
[
ψT(x)K1C1 + ψT(x)K2C2

]
,

ψT(x)C2 = GTψ(x) + +λ2
[
ψT(x)K3C1 + ψT(x)K4C2

]
. (3.6)

Therefore,

C1 ≈ F + λ1 (K1C1 +K2C2) and C2 ≈ G+ λ2 (K3C1 +K4C2) . (3.7)

or,

(I− λ1K1)C1 − λ1K2C2 = F and (I− λ2K3)C2 − λ2K4C1 = G. (3.8)

After replacing ≈ with =, we have a linear system that can be solved with using finite iterative algorithm
for the unknown vectors C1, C2 then by the use of u(x) ≈ C1

Tψ(x), v(x) ≈ C2
Tψ(x) the approximated

solution is given.
The coupled linear system (3.8) can be further written in the form

A1C1 +B1C2 = F and A2C1 +B2C2 = G,

where A1 = I − λ1K1, B1 = −λ1K2, A2 = I − λ2K3, B2 = −λ2K4.

4. Convergence Analysis and Error Estimation

This section provides a brief study of the convergence behavior and error bounds of the proposed
numerical method, ensuring its accuracy and reliability.

4.1. Convergence analysis

In this section, we will discuss the convergence analysis for our proposed numerical approach.

Theorem 3.1 Using our proposed numerical method, Legendre wavelets coupled with a finite iterative
method, the solution of (3.1) converges to ŭ(t) := (ŭ(t), v̆(t)) defined in (3.2).

Proof. Let {ψn,m(t)}n,m be the Legendre wavelets forming an orthonoral basis of L2(R) and let L2(R)

be a Hilbert space. Since ŭ(t) is a vector-valued function with two components in L2(R), we expand it
component -wise

ŭ(t) ≈
2K−1∑
n=1

M−1∑
m=0

c(1)n,mψn,m(t), v̌(t) ≈
2K−1∑
n=1

M−1∑
m=0

c(2)n,mψn,m(t),

where

c(1)n,m =< ŭ(t), ψn,m(t) >, c(2)n,m =< v̌(t), ψn,m(t) > .

Then, the full approximation of the vector function ŭ(t) is



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 6 of 18

ũN (t) ≈
2K−1∑
n=1

M−1∑
m=0

cn,mψn,m(t), where cn,m :=

(
c(1)n,m

c(2)n,m

)
,

and

ψn,m(t) :=

(
ψ
(1)
n,m

ψ
(2)
n,m

)
.

Relabeling (n,m) into a single index j, we write:

ũN (t) =

N∑
j=1

ζjψ (tj) , where ζj =
(
ζ(1)j

ζ(2)j

)
,

here we denote

ζ(1)j :=< ŭ(t), ψj(t) > and ζ(2)j :=< v̌(t), ψj(t) > .

To prove convergence, consider the partial sums:

SN =
N∑
j=1

ζjψj(t).

Let N > M . Then,

‖SN − SM‖2 =

∥∥∥∥∥∥
N∑

j=M+1

ζjψj(t)

∥∥∥∥∥∥
2

=
n∑

j=1

‖ζj‖2 ,

due to the orthonomality of ψj(t).
Hence,

‖SN − SM‖2 = ‖
N∑

j=M+1

(∣∣∣ζ(1)j

∣∣∣2 + ∣∣∣ζ(2)j

∣∣∣2) ,
which converges as M . N → ∞ by Bessel’s inequality. Therefore, {Sn} is a cauchy sequence in

L2(R)× L2(R), and thus converges to some limit ′s(t)′.
To show that s(t) = ŭ(t), consider

< s− ŭ, ψj > = lim
N→∞

< SN − ŭ, ψj >= lim
N→∞

(< SN − ψj > − < ŭ, ψj >)

= ζj − ζj = 0

Thus, the difference s(t)− ŭ(t) is orthogonal to all basis functions ψj(t), and therefore:

s(t) = ŭ(t).

Hence, the approximation series
∑∞

j=1 ζjψj(t) converges to ŭ(t).
Consequently, we have ŭ(t) = s and

∑n
j=1 ζjψ (tj) converges to ŭ(t) in L2(R)×L2(R), completing the proof.



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 7 of 18

4.2. Error Estimation

Suppose that ŭ(t) := (ŭ(t), v̆(t)) is the approximate solution and u(t) := (u(t), v(t)) is the exact
solution, then the error function En(t) is given by the following relation:

En(t) = u(t)− ŭ(t),

hence

ŭ(t) =
2K−1∑
n=1

M−1∑
m=0

cn,mψn,m(t) +Hn(t) = CTψ(t) +Hn(t),

where, Hn(t) is the perturbation term.

Hn(t) = ŭ(t)− CTψ(t) (3.11)

so, its easily to see that

En(t) + CTψ(t) = −Hn(t),

Hence the stability of our proposed method is established through this convergence theorem and error
estimation.

5. Numerical examples

In this section, some numerical examples are provided to illustrate the efficiency and accuracy of our
method, which integrates Legendre wavelets with a proposed finite iterative algorithm. These examples
are drawn from recent existing literature, allowing us to compare the numerical results of our approach
with both exact solutions and those reported in previous studies. All computations were performed using a
program developed in MATLAB R2015a.

Example 1
Consider the system of two linear Fredholm integral equations ([12],[29],[33],[34]):

u1(t) =
t

18
− 17

36
+

∫ 1

x=0

x+ t

3
(u1(x) + u2(x)) dx,

u2(t) = t2 − 19

12
t+ 1 +

∫ 1

t=0
xt (u1(x) + u2(x)) dx, (4.1)

with exact solution (u1(t), u2(t)) =
(
t+ 1, t2 + 1

)
.

This example has been addressed by several researchers. Initially, Babolian et al. [12] solved it using the
Adomian decomposition method. Subsequently, Ramadan et al. [33] tackled it using the triangular basis
functions method with m = 32 triangular basis functions. Najafi et al. [29] also analyzed the problem using
the linear Legendre multi-wavelets method. More recently, Arafa and Ramadan [34] investigated it using
Bernoulli wavelets method with parameters M = 3, k = 2.

Applying our proposed method (LWM), the unknown functions u1(t), u2(t) can be expanded as

u1(t) ≈ CT
1 Ψ(t), u2(t) ≈ CT

2 Ψ(t) (4.2)

where C1, C2 are the unknown 2k−1M, M = 3, k = 2 vectors with
C1 =

[
c1,0 c1,1 c1,2 c2,0 c2,1 c2,2

]T
, C2 =

[
c′1,0 c′1,1 c′1,2 c′2,0 c′2,1c

′
2,2

]T and Ψ(t) is Leg-
endre wavelets for M = 3, k = 2 :



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 8 of 18

Ψ(t) =

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


Ψ1,0 =

√
2

Ψ1,1 =
√
6(4t− 1)

Ψ1,2 =
√
10
(
3
2(4t− 1)2 − 1

2

)
 , 0 ≤ t < 1

2 ,
Ψ2,0 =

√
2

Ψ2,1 =
√
6(4t− 3)

Ψ2,2 =
√
10
(
3
2(4t− 3)2 − 1

2

)
 , 12 ≤ t < 1.

Likewise,

k1(x, t) =
x+ t

3
, k2(x, t) =

x+ t

3
, k3(x, t) = xt, k4(x, t) = xt

f(t) =
t

18
− 17

36
, g(t) = t2 − 19

12
t+ 1,

are also expanded into the LWM as:

k1(x, t) ≈ ΨT (x)K1Ψ(t), k2(x, t) ≈ ΨT (x)K2Ψ(t),

k3(x, t) ≈ ΨT (x)K3Ψ(t), k4(x, t) ≈ ΨT (x)K4Ψ(t),

f(t) ≈ F TΨ(t), g(t) ≈ GTΨ(t). (4.3)

After substituting (4.2), (4.3) into (4.1) we get,

ΨT (t)C1 = F TΨ(t) +
[
ΨT (t)K1C1 +ΨT (t)K2C2

]
,

ΨT (t)C2 = GTΨ(t) +
[
ΨT (t)K3C1 +ΨT (t)K4C2

]
, (4.4)

which can be written in the coupled system of matrix equations,

(I −K1)C1 −K2C2 = F and (I −K3)C2 −K4C1 = G. (4.5)

We can write (4.5) further in the form:

A1C1 +B1C2 = F and A2C1 +B2C2 = G, (4.6)

where,



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 9 of 18

A1 = I −K1 =



0.91667 −0.024056 0 −0.16667 −0.024056 0
−0.024056 1 0 −0.024056 0 0

0 0 1 0 0 0
−0.16667 −0.024056 0 0.75 −0.024056 0
−0.024056 0 0 −0.024056 1 0

0 0 0 0 0 1

 ,

B1 = −K2 =



−0.083333 −0.024056 0 −0.16667 −0.024056 0
−0.024056 0 0 −0.024056 0 0

0 0 0 0 0 0
−0.16667 −0.024056 0 −0.25 −0.024056 0
−0.024056 0 0 −0.024056 0 0

0 0 0 0 0 0

 ,

A2 = −K4 =



−0.03125 −0.018042 0 −0.09375 −0.018042 0
−0.018042 −0.010417 0 −0.054127 −0.010417 0

0 0 0 0 0 0
−0.09375 −0.054127 0 −0.28125 −0.054127 0
−0.018042 −0.010417 0 −0.054127 −0.010417 0

0 0 0 0 0 0

 ,

B2 = I −K3 =



0.96875 −0.018042 0 −0.09375 −0.018042 0
−0.018042 0.98958 0 −0.054127 −0.010417 0

0 0 1 0 0 0
−0.09375 −0.054127 0 0.71875 −0.054127 0
−0.018042 −0.010417 0 −0.054127 0.98958 0

0 0 0 0 0 1

 ,

F =
[
0.34373 0.0056701 0 0.36337 0.0056701 0

]T
,

G =
[
−0.48614 −0.11057 0.013176 0.2799 −0.0085052 0.013176

]T
.

Using our proposed finite iterative algorithm 2.1 to solve the coupled matrix system given in equation
(4.6), we obtain the coefficient vectors

C1 =
[
0.8839, 0.1021, 0, 1.237, 0.1021, 0

]T
and
C2 = [0.766032, 0.051031, 0.0131762, 1.11959, 0.153093, 0.0131762]T .
Substituting these coefficient vectors into a coupled system of equations (4.6) provides the approximate

solutions at the corresponding time points, as shown in Table 1.



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 10 of 18

Table 1: The numerical results for Example 1 with proposed Legendre wavelets for M = 3, k = 2

t Exact solution
(u1(t), u2(t))

Presented Method
(u1(t), u2(t))

Absolute Error

0 1.0 1.0 0.9999999999999529 1.000000000000972 4.7073e− 14 9.7167e− 13

0.1 1.1 1.01 1.099999999999799 1.01000000000094 2.0117e− 13 9.3969e− 13

0.2 1.2 1.04 1.199999999999645 1.040000000000885 3.5505e− 13 8.8507e− 13

0.3 1.3 1.09 1.2999999999999491 1.090000000000808 5.0915e− 13 8.0802e− 13

0.4 1.4 1.16 1.399999999999337 1.160000000000709 6.6303e− 13 7.0877e− 13

0.5 1.5 1.25 1.499999999998751 1.249999999997582 1.2488e− 12 2.4176e− 12

0.6 1.6 1.36 1.599999999998597 1.359999999997396 1.4029e− 12 2.6037e− 12

0.7 1.7 1.49 1.699999999998443 1.489999999997188 1.5568e− 12 2.812e− 12

0.8 1.8 1.64 1.799999999998289 1.639999999996957 1.7109e− 12 3.0429e− 12

0.9 1.9 1.81 1.899999999998135 1.809999999996703 1.8647e− 12 3.2967e− 12

Table 2: Comparison between absolute errors of Example 4.1 for our presented method with M = 3, k = 2, the T.F. method with
m = 32 [33] and Adomian decomposition method proposed in [12].

t TF method [33] Absolute Error Method in [12] Absolute Error Presented Method Absolute Error
Results of u1(t)

0 1.000077 7.74578E − 05 0.988498 1.15020E − 02 0.99999999999995 4.7073e − 14
0.1 1.100097 9.70027E − 05 1.086632 1.33680 E -02 1.0999999999998 2.0117e − 13
0.2 1.200117 1.16548E − 04 1.184766 1.52340E − 02 1.1999999999996 3.5505e − 13
0.3 1.300136 1.36093E − 04 1.282899 1.71010 E -02 1.2999999999995 5.0915e − 13
0.4 1.400156 1.55637E − 04 1.381033 1.89670 E -02 1.3999999999993 6.6303e − 13
0.5 1.500175 1.75182E − 04 1.479167 2.08330 E -02 1.4999999999988 1.2488e − 12
0.6 1.600195 1.94727E − 04 1.577301 2.26990 E -02 1.5999999999986 1.4029e − 12
0.7 1.700214 2.14272E − 04 1.675435 2.45650 E -02 1.6999999999984 1.5568e − 12
0.8 1.800234 2.33817E − 04 1.773569 2.64310 E -02 1.7999999999983 1.7109e − 12
0.9 1.900253 2.53362 E -04 1.9871702 8.71702 E -02 1.8999999999981 1.8647e − 12

Results of u2(t)
t TF method [33] Absolute Error Method in [12] Absolute Error Presented Method Absolute Error
0 1.070363 7.03625E − 02 1.000000 0.00000E + 00 71.000000000001 9.7167e − 13
0.1 1.067708 5.77080E − 02 1.006549 3.45100E − 03 1.0100000000009 9.3969e − 13
0.2 1.086082 4.60824E − 02 1.033099 6.90100E − 03 1.0400000000009 8.8507e − 13
0.3 1.125486 3.54856E − 02 1.079648 1.03520E − 02 1.0900000000008 8.0802e − 13
0.4 1.185918 2.59177E − 02 1.146198 1.38020E − 02 1.1600000000007 7.0877e − 13
0.5 1.267379 1.73787E − 02 1.232747 1.72530E − 02 1.2499999999976 2.4176e − 12
0.6 1.369869 9.86860E − 03 1.339296 2.07040E − 02 1.35999999999974 2.6037e − 12
0.7 1.493387 3.38735E − 03 1.465846 2.41540E − 02 1.4899999999972 2.812e − 12
0.8 1.637935 2.06503E − 03 1.612695 2.73050E − 02 1.639999999997 3.0429e − 12
0.9 1.803511 6.48854E − 03 1.778945 3.10550E − 02 1.8099999999967 3.2967e − 12

As shown in Tables 2 and 3, our proposed method yields more accurate results compared to those
reported in references [12], [33], and [29]. Furthermore, Table 4 indicates that the accuracy of our method is
nearly equivalent to that of the approach presented in reference [34]. In addition to its superior or compara-
ble accuracy, our method is also highly efficient, as it eliminates the need for the computationally intensive
matrix inversion when determining the coefficient vectors C1, C2.

Example 2
Consider the system of two linear Fredholm integral equations ([12],[25,[33]).

u1(t) = t− 5
18 + 1

3

∫ 1
x=0 (u1(x) + u2(x)) dx,

u2(t) = t2 − 5
12 + 1

2

∫ 1
x=0 (u1(x) + u2(x)) dx,

(4.7)

Table 3: Comparison between absolute errors for results of Example 1 for our presented method with M = 3, k = 2 and presented
method [29].

t Method in [29] Absolute Error Presented Method Absolute error Method in [29] Absolute error Presented Method Absolute error
Results of u1(t) Results of u2(t)

0 0 4.7073e − 14 8.437524563517452E − 01 9.7167e − 13
0.1 0 2.0117e − 13 4.583333669999856E − 03 9.3969e − 13
0.2 0 3.5505e − 13 4.166662600000315E − 04 8.8507e − 13
0.3 0 5.0915e − 13 4.166661900000257E − 04 8.0802e − 13
0.4 0 6.6303e − 13 4.583333879999874E − 03 7.0877e − 13
0.5 0 1.2488e − 12 1.041666597500002E − 02 2.4176e − 12
0.6 0 1.4029e − 12 4.583334199999900E − 03 2.6037e − 12
0.7 0 1.5568e − 12 4.166656999999852E − 04 2.812e − 12
0.8 0 1.7109e − 12 4.166658000002155E − 04 3.0429e − 12
0.9 0 1.8647e − 12 4.583334299999686E − 03 3.2967e − 12



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 11 of 18

Table 4: Comparison between absolute errors for results of Example 1 for our presented method with M = 3, k = 2 and presented
method in [34].

t Method in[34] Absolute error Presented Method Absolute error Method in[34] Absolute error Presented Method Absolute error
Results of u1(t) Results of u2(t)

0 4.7073e − 14 4.7073e − 14 9.7167e − 13 9.7167e − 13
0.1 2.0117e − 13 2.0117e − 13 9.3969e − 13 9.3969e − 13
0.2 3.5505e − 13 3.5505e − 13 8.8507e − 13 8.8507e − 13
0.3 5.0915e − 13 5.0915e − 13 8.0802e − 13 8.0802e − 13
0.4 6.6303e − 13 6.6303e − 13 7.0877e − 13 7.0877e − 13
0.5 1.2488e − 12 1.2488e − 12 2.4176e − 12 2.4176e − 12
0.6 1.4029e − 12 1.4029e − 12 2.6037e − 12 2.6037e − 12
0.7 1.5568e − 12 1.5568e − 12 2.812e − 12 2.812e − 12
0.8 1.7109e − 12 1.7109e − 12 3.0429e − 12 3.0429e − 12
0.9 1.8647e − 12 1.8647e − 12 3.2967e − 12 3.2967e − 12

with exact solution (u1(t), u2(t)) =
(
t, t2

)
.

This example has been addressed by several researchers. Initially, Ramadan et al. [33] tackled it using
the Triangular Basis Functions Method with m = 32 and m = 10 using triangular basis
functions, also authors in [25] analyzed this problem using Triangular functions method, besides , E. Babolian
et al. [12] investigated the same problem using Adomian decomposition method.

By taking M = 3, k = 2, applying our proposed method (LWM), the unknown functions u1(t), u2(t)
can be expanded as

u1(t) ≈ CT
1 Ψ(t), u2(t) ≈ CT

2 Ψ(t), (4.8)

where C1, C2 are the unknown 2k−1M,M = 3, k = 2 vectors with
C1 =

[
c1,0 c1,1 c1,2 c2,0 c2,1 c2,2

]T , C2 =
[
c′1,0 c′1,1 c′1,2 c′2,0 c′2,1 c′2,2

]T and Ψ(t) is Legen-
dre wavelets for M = 3, k = 2 :

Ψ(t) =

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


Ψ1,0 =

√
2

Ψ1,1 =
√
6(4t− 1)

Ψ1,2 =
√
10
(
3
2(4t− 1)2 − 1

2

)
 , 0 ≤ t < 1

2 ,
Ψ2,0 =

√
2

Ψ2,1 =
√
6(4t− 3)

Ψ2,2 =
√
10
(
3
2(4t− 3)2 − 1

2

)
 , 12 ≤ t < 1.

Likewise,

k1(x, t) =
1

3
, k2(x, t) =

1

3
, k3(x, t) =

1

2
, k4(x, t) =

1

2
, f(t) = t− 5

18

g(t) = t2 − 5
12 , are also expanded into the LWM as:

k1(x, t) ≈ ΨT (x)K1Ψ(t), k2(x, t) ≈ ΨT (x)K2Ψ(t),

k3(x, t) ≈ ΨT (x)K3Ψ(t), k4(x, t) ≈ ΨT (x)K4Ψ(t),

f(t) ≈ F TΨ(t), g(t) ≈ GTΨ(t). (4.9)

After substituting (4.8), (4.9) into (4.7) we get,

ΨT (t)C1 = F TΨ(t) +
[
ΨT (t)K1C1 +ΨT (t)K2C2

]
,

ΨT (t)C2 = GTΨ(t) +
[
ΨT (t)K3C1 +ΨT (t)K4C2

]
, (4.10)

which can be written in the coupled system of matrix equations,

(I −K1)C1 −K2C2 = F and (I −K3)C2 −K4C1 = G. (4.11)



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 12 of 18

We can write (4.11) further in the form:

A1C1 +B1C2 = F and A2C1 +B2C2 = G, (4.12)

where,

A1 = I −K1 =



0.83333 0 0 −0.16667 0 0
0 1 0 0 0 0
0 0 1 0 0 0

−0.16667 0 0 0.83333 0 0
0 0 0 0 1 0
0 0 0 0 0 1

 ,

B1 = −K2 =



−0.16667 0 0 −0.16667 0 0
0 0 0 0 0 0
0 0 0 0 0 0

−0.16667 0 0 −0.16667 0 0
0 0 0 0 0 0
0 0 0 0 0 0

 ,

A2 = −K4 =



−0.25 0 0 −0.25 0 0
0 0 0 0 0 0
0 0 0 0 0 0

−0.25 0 0 −0.25 0 0
0 0 0 0 0 0
0 0 0 0 0 0

 ,

B2 = I −K3 =



0.75 0 0 −0.25 0 0
0 1 0 0 0 0
0 0 1 0 0 0

−0.25 0 0 0.75 0 0
0 0 0 0 1 0
0 0 0 0 0 1

 ,

F =
[
−
√
2

72 ,
√
6

24 , 0,
17

√
2

72 ,
√
6

24 , 0
]T

,

G =
[
−
√
2

6 ,
√
6

48 ,
√
10

240 ,
√
2

12 ,
√
6

16 ,
√
10

240

]T
.

Using our proposed finite iterative algorithm 2.1 to solve the coupled matrix system given in equation
(4.12), we obtain the coefficient vectors

C1 = [0.17679, 0.10206, 0, 0.53034, 0.10206, 0]T ,

and

C2 = [0.058937, 0.051031, 0.013176, 0.41249, 0.15309, 0.013176]T .

Substituting these coefficient vectors into coupled system of equations (4.12) provides the approximate
solutions at the corresponding time points, as shown in Table 5 below.

As shown in Tables 6 and 7, our proposed method yields more accurate results compared to those
reported in references [33], [12], and [25]. In addition to its superior or comparable accuracy, our method
is also highly efficient, as it eliminates the need for the computationally intensive matrix inversion when
determining the coefficient vectors C1, C2

Example 3
Consider the system of two coupled linear Fredholm integral equations [30]



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 13 of 18

Table 5: The numerical results for Example 2 with proposed Legendre wavelets for M = 3, k = 2

t Exact solution Presented Method (u1(t), u2(t)) Absolute Error
0 0 0 0.00002389254361656518 0.00001576411199376885 2.3893e− 05 1.5764e− 05

0.1 0.1 0.01 0.1000218618029359 0.01001620490521751 2.1862e− 05 1.6205e− 05

0.2 0.2 0.04 0.2000198310622553 0.04001640751462307 1.9831e− 05 1.6408e− 05

0.3 0.3 0.09 0.3000178003215747 0.09001637194021046 1.78e− 05 1.6372e− 05

0.4 0.4 0.16 0.400015769580894 0.1600160981819797 1.577e− 05 1.6098e− 05

0.5 0.5 0.25 0.5000190975197177 0.2500227380709804 1.9098e− 05 2.2738e− 05

0.6 0.6 0.36 0.6000170667790372 0.3600201683276863 1.7067e− 05 2.0168e− 05

0.7 0.7 0.49 0.7000150360383566 0.490017360400574 1.5036e− 05 1.736e− 05

0.8 0.8 0.64 0.8000130052976759 0.6400143142896436 1.3005e− 05 1.4314e− 05

0.9 0.9 0.81 0.9000109745569953 0.8100110299948949 1.0975e− 05 1.103e− 05

Table 6: Comparison between absolute errors of Example 4.2 for our presented method with M = 3, k = 2, the T.F. method with
m = 32 [33] and when m = 10 [33]

t TF method [33] (m = 10) Absolute Error TF method [33] ( m = 32 ) Absolute Error Presented Method Absolute Error
Results of u1(t)

0 0.00333333 2.49088e − 003 0.00032552 3.25520839e − 04 0.00002389254361656518 2.3893e − 05
0.1 0.10333333 2.49083e − 003 0.10032552 3.25520839e − 04 0.1000218618029359 2.1862e − 05
0.2 0.20333333 2.49078e − 003 0.20032552 3.25520839e − 04 0.2000198310622553 1.9831e − 05
0.3 0.30333333 2.49073e − 003 0.30032552 3.25520839e − 04 0.3000178003215747 1.78e − 05
0.4 0.40333333 2.49068e − 003 0.40032552 3.25520839e − 04 0.400015769580894 1.577e − 05
0.5 0.50333333 2.49063e − 003 0.50032552 3.25520839e − 04 0.5000190975197177 1.9098e − 05
0.6 0.60333333 2.49058e − 003 0.60032552 3.25520839e − 04 0.6000170667790372 1.7067e − 05
0.7 0.70333333 2.49053e − 003 0.70032552 3.25520839e − 04 0.7000150360383566 1.5036e − 05
0.8 0.80333333 2.49047e − 003 0.80032552 3.25520839e − 04 0.8000130052976759 1.3005e − 05
0.9 0.90333333 2.49042e − 003 0.90032552 3.25520839e − 04 0.9000109745569953 1.0975e − 05

Results of u2(t)
t TF method [33] For m=10 Absolute Error TF method [33] Absolute Error Presented Method Absolute
0 0.00500000 4.15743e − 003 -0.00146484 1.46484374e − 03 0.00001576411199376885 1.5764e − 05
0.1 0.01500000 4.15743e − 003 0.00908203 9.17968744e − 04 0.01001620490521751 1.6205e − 05
0.2 0.04500000 4.15741e − 003 0.03955078 4.49218744e − 04 0.04001640751462307 1.6408e − 05
0.3 0.09500000 4.15739e − 003 0.08994141 5.85937439e − 05 0.09001637194021046 1.6372e − 05
0.4 0.16500000 4.15735e − 003 0.16025391 2.53906256e − 04 0.1600160981819797 1.6098e − 05
0.5 0.25500000 4.15731e − 003 0.25048828 4.88281256e − 04 0.2500227380709804 2.2738e − 05
0.6 0.36500000 4.15725e − 003 0.36064453 6.44531256e − 04 0.3600201683276863 2.0168e − 05
0.7 0.49500000 4.15718e − 003 0.49072266 7.22656256e − 04 0.490017360400574 1.736e − 05
0.8 0.64500000 4.15711e − 003 0.64072266 7.22656256e − 04 0.6400143142896436 1.4314e − 05
0.9 0.81500000 4.15702e − 003 0.81064453 6.44531256e − 04 0.8100110299948949 1.103e − 05

u1(t) = sin(t)− cos(1) + sin(1)− t sin(1) +

∫ 1

x=0
[(t− s)u1(s) + tsu2(s)] ds,

u2(t) = cos(t)− (1− cos(1))t2 + cos(1)− 3 sin(1)− t sin(1) + 1+∫ 1

x=0

[(
t2 + 2s

)
u1(s) + (s+ t)u2(s)

)]
ds, (4.13)

with exact solution (u1(t), u2(t)) = (sin(t), cos(t)).
Noting that, Zaffer Elahi et al. [30] tackled this system using Laguerre method.

By taking M = 3, k = 2, applying our proposed method (LWM), the unknown functions u1(t), u2(t)
can be expanded as

u1(t) ≈ CT
1 Ψ(t), u2(t) ≈ CT

2 Ψ(t) (4.14)

where C1, C2 are the unknown 2k−1M,M = 3, k = 2 vectors with
C1 =

[
c1,0 c1,1 c1,2 c2,0 c2,1 c2,2

]T and C2 =
[
c′1,0 c′1,1 c′1,2 c′2,0 c′2,1 c′2,2

]T and Ψ(t) is
Legendre wavelets for M = 3, k = 2 :



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 14 of 18

Table 7: Comparison between absolute errors for results of Example 4.2 for our presented method with M = 3, k = 2 and Absolute
Error for TF method [25] and Absolute Error for Adomian method in [12]

t AE for TF method [25] AE in [12] AE for our presented AE in [12] AE for our presented
Results of u1(t) Results of u2(t)

0 2.170e− 04 2.30e− 02 2.3893e− 05 4.43e− 02 1.5764e− 05

0.1 2.170e− 04 2.30e− 02 2.1862e− 05 4.43e− 02 1.6205e− 05

0.2 2.170e− 04 2.30e− 02 1.9831e− 05 4.43e− 02 1.6408e− 05

0.3 2.170e− 04 2.30e− 02 1.78e− 05 4.43e− 02 1.6372e− 05

0.4 2.170e− 04 2.30e− 02 1.577e− 05 4.43e− 02 1.6098e− 05

0.5 2.170e− 04 2.30e− 02 1.9098e− 05 4.43e− 02 2.2738e− 05

0.6 2.170e− 04 2.30e− 02 1.7067e− 05 4.43e− 02 2.0168e− 05

0.7 2.170e− 04 2.30e− 02 1.5036e− 05 4.43e− 02 1.736e− 05

0.8 2.170e− 04 2.30e− 02 1.3005e− 05 4.43e− 02 1.4314e− 05

0.9 2.170e− 04 2.30e− 02 1.0975e− 05 4.43e− 02 1.103e− 05

Ψ(t) =

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


Ψ1,0 =

√
2

Ψ1,1 =
√
6(4t− 1)

Ψ1,2 =
√
10
(
3
2(4t− 1)2 − 1

2

)
 , 0 ≤ t < 1

2 ,
Ψ2,0 =

√
2

Ψ2,1 =
√
6(4t− 3)

Ψ2,2 =
√
10
(
3
2(4t− 3)2 − 1

2

)
 , 12 ≤ t < 1.

Likewise,

k1(s, t) = (t− s), k2(s, t) = ts, k3(s, t) =
(
t2 + 2s

)
k4(s, t) = (s+ t), f(t) = sin(t)− cos(1) + sin(1)− t sin(1)

g(t) = cos(t)− (1− cos(1))t2 + cos(1)− 3 sin(1)− t sin(1) + 1,

are also expanded into the LWM as:

k1(s, t) ≈ ΨT (s)K1Ψ(t), k2(s, t) ≈ ΨT (s)K2Ψ(t)

k3(s, t) ≈ ΨT (s)K3Ψ(t), k4(s, t) ≈ ΨT (s)K4Ψ(t)

f(t) ≈ F TΨ(t), g(t) ≈ GTΨ(t) (4.15)

After substituting (4.14), (4.15) into (4.13) we get,

ΨT (t)C1 = F TΨ(t) +
[
ΨT (t)K1C1 +ΨT (t)K2C2

]
ΨT (t)C2 = GTΨ(t) +

[
ΨT (t)K3C1 +ΨT (t)K4C2

]
(4.16)

which can be written in the coupled system of matrix equations,

(I −K1)C1 −K2C2 = F and (I −K3)C2 −K4C1 = G. (4.17)

We can write (4.3.5) further in the form,

A1C1 +B1C2 = F and A2C1 +B2C2 = G, (4.18)



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 15 of 18

where,

A1 = I −K1 =



1 0.072169 0 0.25 0.0721690 0
−0.072169 1 0 −0.072169 0 0

0 0 1 0 0 0
−0.25 0.072169 0 1 0.072169 0

−0.072169 0 0 −0.072169 1 0
0 0 0 0 0 1

 ,

B1 = −K2 =



−0.03125 −0.018042 0 −0.09375 −0.018042 0
−0.018042 −0.010417 0 −0.054127 −0.010417 0

0 0 0 0 0 0
−0.09375 −0.054127 0 −0.28125 −0.054127 0
−0.018042 −0.010417 0 −0.054127 −0.010417 0

0 0 0 0 0 0

 ,

A2 = −K4 =



−0.29167 −0.14434 0 −0.79167 −0.14434 0
−0.036084 0 0 −0.036084 0 0
−0.0093169 0 0 −0.0093169 0 0
−0.54167 −0.14434 0 −1.0417 −0.14434 0
−0.10825 0 0 −0.10825 0 0
−0.0093169 0 0 −0.0093169 0 0

 ,

B2 = I −K3 =



0.75 −0.072169 0 −0.5 −0.072169 0
−0.072169 1 0 −0.072169 0 0

0 0 1 0 0 0
−0.5 −0.072169 0 0.25 −0.072169 0

−0.072169 0 0 −0.072169 1 0
0 0 0 0 0 1

 ,

F =
[
0.23733 0.01239 −0.0016227 0.24369 −0.01167 −0.0044707

]T
,

G =
[
−0.1937 −0.13443 −0.012412 −0.81973 −0.22539 −0.010856

]T
.

Using our proposed finite iterative algorithm 2.1 to solve the coupled matrix system given in equation
(4.18), we obtain the coefficient vectors

C1 = [0.1731, 0.0983,−0.0016, 0.477, 0.0742,−0.0045]T ,

and

C2 = [0.678,−0.0251,−0.0064, 0.512,−0.0691,−0.0048]T .

Substituting these coefficient vectors into coupled system of equations (4.18) provides the approximate
solutions at the corresponding time points, as shown in Table 8

Table 8: The numerical results for Example 3 with proposed Legendre wavelets for M = 3, k = 2

T Exact solution (u1(t), u2(t)) Presented Method (u1(t), u2(t)) Absolute Error (u1(t), u2(t))
0 0 1.0 -0.001044118325073061 1.000080410807739 0.00104 8.04e-5
0.2 0.19866933 0.98006658 0.198869642776424 0.9800382076887642 2.0e− 4 2.84e− 5
0.4 0.38941834 0.92106099 0.3890688869058838 0.9211379366816407 3.49e− 4 7.69e− 5
0.6 0.56464247 0.82533561 0.5649593759244289 0.8250260313600629 3.17e− 4 3.1e− 4
0.8 0.71735609 0.69670671 0.7171916068020021 0.6969041261080368 1.64e− 4 1.97e− 4
1 0.84147098 0.54030231 0.8421017586957205 0.5396386699398988 6.31e− 4 6.64e− 4



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 16 of 18

Table 9: Comparison between absolute errors of Example 3 for our presented method with M = 3, k = 2, and Laguerre method in
[30] for N = 2

t AE in [30] taking N = 2 AE for Presented Method AE in [30] AE for Presented Method
Results of u1(t) Results of u2(t)

0 6.12984 e -4 0.00104 1.71016 e -3 8.04e− 5

0.2 7.9168 e -3 2.0e− 4 5.33819 e -3 2.84e− 5

0.4 4.36011 e -3 3.49e− 4 4.18645 e -3 7.69e− 5

0.6 2.45251 e -3 3.17e− 4 6.07311 e -4 3.1e− 4

0.8 5.53543 e -3 1.64e− 4 1.58297 e -3 1.97e− 4

1 1.19956 e -3 6.31e− 4 2.74364 e -3 6.64e− 4

As shown in Table 9, our proposed method yields more accurate results compared to obtained results in
[30]. In addition to its superior or comparable accuracy, our method is also highly efficient, as it eliminates
the need for the computationally intensive matrix inversion when determining the coefficient vectors C1, C2.

6. Conclusion

This study presents an efficient and direct iterative method for addressing one-dimensional Fredholm
integral equations of the second order, utilizing Legendre wavelet functions and a finite iterative framework.
The method transforms integral equations into systems of algebraic matrix equations, providing a practical
and efficient approach for calculating approximation solutions. The method circumvents the necessity of
inverting block matrices, thereby improving its precision and computational efficiency. Numerous numerical
cases exemplify the method’s performance, showcasing its potential for diverse applications in scientific and
technical domains.

Dedication

The first author, Mohamed A. Ramadan, dedicates this work to Professor Mohamed Asaad, Profes-
sor Emeritus at Cairo University, on his 80th birthday. His pioneering research in abstract algebra and
finite groups is a true inspiration. Though my research lies outside his field, his generous support and
encouragement have been invaluable.

References

[1] G. B. Arfken and H. J. Weber. Mathematical Methods for Physicists. Elsevier Academic Press, 6th
edition, 2005.

[2] L. Debnath and D. Bhatta. Integral Transforms and Their Applications. Chapman and Hall/CRC, 3rd
edition, 2014.

[3] A. D. Polyanin and A. V. Manzhirov. Handbook of Integral Equations. Chapman and Hall/CRC, 2nd
edition, 2008.

[4] Z. Elahi, G. Akram, and S. S. Siddiqi. Numerical solutions for solving special eighth-order linear
boundary value problems using legendre galerkin method. Mathematical Sciences, 10(4):201–209, 2016.

[5] Z. Elahi, G. Akram, and S. S. Siddiqi. Numerical solutions for solving special tenth-order linear bound-
ary value problems using legendre galerkin method. Mathematical Sciences Letters, 7(1):27–35, 2018.

[6] N. Negarchi and K. Nouri. Numerical solution of volterra-fredholm integral equations using the collo-
cation method based on a special form of the müntz-legendre polynomials. Journal of Computational
and Applied Mathematics, 344:15–24, 2018.

[7] S. Nemati. Numerical solution of volterra-fredholm integral equations using legendre collocation method.
Journal of Computational and Applied Mathematics, 278:29–36, 2015.



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 17 of 18

[8] B. Yilmaz and Y. Cetin. Numerical solutions of the fredholm integral equations of second type. NTM-
SCI, 3(5):284–292, 2017.

[9] Mohamed A. Ramadan and Mohamed R. Ali. Solution of integral and integro-differential equations
system using hybrid orthonormal bernstein and block-pulse functions. Journal of Abstract and Compu-
tational Mathematics, NTMSCI, 2(1):35–48, 2017.

[10] Mohamed A. Ramadan and Mohamed R. Ali. An efficient hybrid method for solving fredholm integral
equations using triangular functions. New Trends in Mathematical Sciences, 5(1):213–224, 2017.

[11] M. Ramadan, K. Raslan, A. Hadhoud, and M. Nassar. Rational chebyshev-based schemes for handling
high-order linear integro-differential equations. 2016. Unpublished manuscript description.

[12] E. Babolian, J. Biazar, and A. R. Vahidi. The decomposition method applied to systems of fredholm
integral equations of the second kind. Applied Mathematics and Computation, 148(2):443–452, 2004.

[13] E. Babolian, Z. Masouri, and S. Hatamzadeh-Varmazyar. A direct method for numerically solving
integral equations system using orthogonal triangular functions. International Journal of Industrial
Mathematics, 1(2):135–145, 2009.

[14] A. Jafarian and S. Measoomy Nia. Utilizing feed-back neural network approach for solving linear
fredholm integral equations system. Applied Mathematical Modelling, 37(7):5027–5038, 2013.

[15] A. Jafarian, S. A. Measoomy Nia, A. K. Gulmankhaneh, and D. Baleanu. Numerical solution of linear
integral equations system by using the bernstein collocation method. Advances in Difference Equations,
2013:123, 2013.

[16] Z. Mahmoodi. Collocation method for solving systems of fredholm and volterra integral equations.
International Journal of Computer Mathematics, 91(8):1802–1816, 2013.

[17] M. Alipour, D. Baleanu, and K. Karimi. Spectral method based on bernstein polynomial for coupled
system of fredholm integral equations. Applied and Computational Mathematics, 15(2):212–219, 2016.

[18] Y. Huang, M. Fang, and X.-F. Li. Approximate solution of a system of linear integral equations by the
taylor expansion method. International Journal of Computer Mathematics, 86(5):924–937, 2009.

[19] Mohamed A. Ramadan, Mokhtar A. Abdel-Naby, and Ahmed M. Bayoumi. Iterative algorithm for
solving a class of general sylvester-conjugate matrix equation. Journal of Applied Mathematics and
Computing, 44:99–118, 2014.

[20] Mohamed A. Ramadan and Talaat S. El-Danaf. Solving the generalized coupled sylvester matrix
equations over generalized bisymmetric matrices. Transactions of the Institute of Measurement and
Control, 37(3):291–316, 2015.

[21] K. Maleknejad, H. Safdari, and M. Nouri. Numerical solution of an integral equations system of the
first kind by using an operational matrix with block pulse functions. International Journal of Systems
Science, 42(1):195–199, 2011.

[22] M. A. Ramadan, H. S. Osheba, and A. R. Hadhoud. A highly efficient and accurate finite iterative
method for solving linear two-dimensional fredholm fuzzy integral equations of the second kind using
triangular functions. Mathematical Problems in Engineering, 2020:Article ID 2028763, 1–16, 2020.

[23] K. Maleknejad, M. Shahrezaee, and H. Khatami. Numerical solution of integral equations system of
the second kind by block-pulse functions. Applied Mathematics and Computation, 166:15–24, 2005.

[24] H. Almasieh and M. Roodaki. Triangular functions method for the solution of fredholm integral equa-
tions system. Ain Shams Engineering Journal, 3(4):411–416, 2012.

[25] A. Golbabai and B. Keramati. Easy computational approach to solution of system of linear fredholm
integral equations. Chaos, Solitons and Fractals, 38(2):568–574, 2008.

[26] P. K. Sahu and S. Saha Ray. Numerical solutions for the system of fredholm integral equations of
second kind by a new approach involving semiorthogonal b-spline wavelet collocation method. Applied
Mathematics and Computation, 234:368–379, 2014.

[27] H. S. Najafi, H. Aminikhah, and S. A. Edalatpanah. Linear legendre multiwavelets methods for solving
systems of fredholm integral equations. Mathematical Reports, 18:41–50, 2016.

[28] Z. Elahi, S. S. Siddiqi, and G. Akram. Laguerre method for solving linear system of fredholm integral
equations. International Journal of Computer Mathematics, 98(11):2175–2185, 2021.



M. A. Ramadan et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6725 18 of 18

[29] Z. Elahi, G. Akram, and S. S. Siddiqi. Laguerre approach for solving system of linear fredholm integro-
differential equations. Mathematical Sciences, 12(3):185–195, 2018.

[30] Z. Elahi, G. Akram, and S. S. Siddiqi. Use of bessel polynomials for solving differential difference
equations. Arab Journal of Basic and Applied Sciences, 26(1):23–29, 2019.

[31] M. Ramadan, H. Oshaba, and R. Kharabsheh. Triangular functions-based method for the solution
of system of linear fredholm integral equations via an efficient finite iterative algorithm. Journal of
Intelligent and Fuzzy Systems, 38(3):2847–2858, 2020.

[32] H. M. Arafa and M. A. Ramadan. Bernoulli wavelet method for numerical solution of linear system of
fredholm integral equation of the second kind. Alexandria Engineering Journal, 77:63–74, 2023.

[33] M.A. Ramadan H.M. Arafa. Bernoulli wavelet method for numerical solution of linear system of fred-
holm integral equation of the second kind. Alexandria Engineering Journal, 77:64–75, 2023.