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

Derivative-Free Method for Nonlinear Systems with
Some Real-life Applications

Yanal Al-Shorman1,∗, Abdullah Alsoboh2, Obadah Said Solaiman1,
Ishak Hashim1,3

1 School of Mathematical Sciences, Faculty of Science & Technology, Universiti Kebangsaan
Malaysia, 43600 Bangi, Selangor, Malaysia
2 Department of Basic and Applied Sciences, College of Applied and Health Sciences,
A’Sharqiyah University, Post Box No. 42, Post Code No. 400, Ibra, Sultanate of Oman
3 Nonlinear Dynamics Research Center (NDRC), Ajman University, Ajman,
P.O. Box 346, United Arab Emirates

Abstract. This study implements a very effective derivative-free approach for root determina-
tion in the context of nonlinear equations. The use of Homeier’s third-order technique, which
is based on the application of Newton’s theorem to the inverse function, as well as the develop-
ment of a new category of Newton-type procedures that exhibit cubic convergence, is the basis
of this study. A variety of non-linear problems, such as non-linear equations and non-linear
systems of equations, have been examined to examine the efficacy of the proposed method in
contrast to other widely used derivative-free approaches. The results suggest that the proposed
method outperforms the currently examined techniques published. The scheme given in this
study demonstrates a reduced number of iterations acquired in attaining the intended result,
exhibiting a convergence order of around 2.5. This convergence order surpasses the convergence
order observed in comparable approaches. Moreover, by following the prescribed procedure of
the Broyden technique while implementing the recommended approach in the solution of nonlin-
ear equation systems, the problem pertaining to the Jacobian can be circumvented. Hence, the
technique suggested may be regarded as a prominent approach that achieves rapid convergence in
determining the roots of non-linear equations without the use of derivatives. In practice, we em-
ploy the novel approach to address practical challenges such as the mixed Hammerstein integral
equation, the discretized Poisson equation, and the Weierstrass method. Through comparisons
and illustrative examples, it is evident that the proposed method proves to be effective and on
par with existing techniques of equivalent order.

2020 Mathematics Subject Classifications: 65H10, 90C56

Key Words and Phrases: Homeier method, Secant method, Nonlinear equations, Derivative-
free methods, Broyden’s method, Iterative methods, Order of convergence, Hammerstein integral
equation

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6262

Email addresses: P114391@siswa.ukm.edu.my (Y. Al-Shorman),
abdullah.alsoboh@asu.edu.om (A. Alsoboh)

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



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 2 of 20

1. Introduction

Solving systems of nonlinear equations has been a central focus in numerical anal-
ysis, leading to the development of numerous iterative techniques [1–4] and references
therein. Among these, Newton’s method remains one of the most widely used approaches
for finding a single root of a nonlinear equation [5–19]. An enhancement of this method
was proposed by Homeier [20, 21], which achieves cubic convergence, offering improved
performance over the classical Newton method.

Both Newton’s and Homeier’s methods require the evaluation of the derivative of the
function, which can limit their applicability in scenarios where the derivative is difficult
or impossible to compute. This issue commonly arises in practical modeling situations, as
discussed in [22]. To address such limitations, optimization-based approaches have been
proposed as alternative strategies for solving nonlinear equations and systems. Over the
past decade, a variety of innovative optimization techniques have emerged [23].

Additionally, several derivative-free iterative schemes have been introduced for solv-
ing systems of nonlinear equations, offering viable alternatives when derivative informa-
tion is unavailable [24, 25]. One such method is the secant approach, which approximates
the derivative in Newton’s method using finite differences. This technique achieves a
convergence order of approximately 1.618 and requires two function evaluations per it-
eration.

However, the existing pace appears to be insufficient in light of the contemporary
global landscape, which requires more efficient algorithms to meet the escalating de-
mands driven by the rapid advancements in technology. Consequently, we were moti-
vated to explore an approach that addresses the present requirements by substituting
the derivatives in the fast-converging Homeier method with suitable difference approxi-
mations. This is necessitated by the absence of robust numerical techniques that do not
rely on derivatives for approximating the solutions of a non-linear equation. This study
presents a methodology for developing an optimal algorithm that can effectively solve
nonlinear equations numerically without relying on derivatives. The proposed approach
aims to optimize both the convergence rate and the number of function evaluations to
achieve high efficiency. Subsequently, an account is provided of the verification process
employed to ascertain the suitability of the approach for solving non-linear equations.
Ultimately, the text elucidates the manner in which the said approach was modified to
include non-linear equation systems, yielding commendable outcomes.

In this paper, we introduce an iterative approach to addressing nonlinear models.
The novel approach represents an adaptation of the Homeier method. The utilization
of the composition methodology, including both forward and backward approximations,
is employed in the development of our novel method. The assignment of tasks in this
investigation is outlined below. The early findings are presented in Section 2. Section
3 illustrates the revised system. The analysis of the convergence sequence of the new



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 3 of 20

scheme is addressed in Section 4. Section 5 of the document provides a detailed explana-
tion of numerical examples and comparisons. Section 6 presents the real-life application.
Finally, the conclusion is presented in Section 7.

2. Preliminary Results

Definition 1. Given that the iterative sequence xn : n = 1, 2, . . . converges to x∗ to
be formed utilizing a numerical scheme. Suppose that there exists a constant c ≥ 0.
Therefore, integers n0 ≥ 0 and ρ ≥ 0 in which for all n > n0, the inequality provided
below holds for any vector norm ∥.∥.

∥xn+1 − x∗ ∥ ≤ c ∥xn − x∗∥ρ (1)

Hence, it is well established that the iterative scheme converges to the solution x∗ having
a convergence order of ρth.

Definition 2. Let x∗ be a root of equation F (x) = 0. Let xn−1, xn as well as xn+1

be consecutive iterations closer to the root x∗ formed by an iterative scheme. Therefore,
the Computational Order of Convergence (COC) ρ with regards to the numerical
algorithm or the iterative scheme may be approximated by

ρ ≈
ln ∥xn+1 − x∗ ∥

∥xn − x∗ ∥

ln ∥xn − x∗ ∥
∥xn−1 − x∗ ∥

(2)

2.1. The Homeier Method in solving systems of nonlinear equations

Provided that F denotes a vector-valued function with non-zero derivatives defined on
the set D⊂ Rn and X, X0∈D, the application of Homeier method’s extension in solving
nonlinear equation’s systems can be articulated as follows

Xn+1=Xn−
1

2

[
F ′ (Xn )

−1
+F ′ (X∗

n+1 )
−1

]
(F (Xn )) (3)

where
X∗

n+1=Xn−F ′ (Xn )
−1

(F (Xn )) (4)

Here, n = 0, 1, 2, · · · and Xn is the nth iterate.

2.2. Secant Method for Systems of Nonlinear Equations

A particular issue arises when using the secant approach to solve systems of nonlinear
equations, specifically in relation to the Jacobian. Here, the secant technique for sys-
tems involves the inclusion of a vector in the denominator. It is customary to employ
the widely used secant approximation method that Broyden [26] proposed to address
the difficulty of computing the inverse of a vector. The algorithm refers to Newton’s ap-
proach. Nevertheless, the subsequent approximation replaces the analytic Jacobian. The
procedure in question is referred to as Broyden’s method, as described by Heenatigala
et al. [27].



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 4 of 20

2.3. Broyden’s Method for Systems of Nonlinear Equations

2.3.1. Application of DFH for systems of nonlinear equations

An initial estimate given by X0 ∈ Rn is selected along with a non-singular initial
matrix A0 ∈ Rn×n. By setting k := 0, we repeat the sequence of steps given below until
∥F (Xk)∥ < tolerance

(i) Solve AkIk = F (X∗
k+1)− F (Xk) for Ik & X∗

k+1 from Broyden’s method.

(ii) Xk+1 = Xk + Ik

(iii) Yk = F (Xk+1)− F (Xk)

(iv) Ak+1 = Ak +
Yk−AkIk

ItkIk
Itk

3. Derivation of Derivative-Free Homeier (DFH) from Homeier
method

In this section, we will explain how our new method was derived. In this context,
a combination of backward and forward difference approximations is employed in a
judicious manner within a single formula to provide a satisfactory outcome. Consider
the Homeier Method

xn+1 = xn − f (xn)

2

[
1

f ′ (xn)
+

1

f ′ (yn)

]
(5)

in which yn = xn − f(xn)
f ′(xn)

with n = 0, 1, 2, · · · .

Replacing the forward difference approximation for

f ′ (xn) ≈
f (xn+1)− f (xn)

xn+1 − xn
(6)

as well as the backward difference approximation for

f ′ (xn+1) ≈
f (xn+1)− f (xn)

xn+1 − xn
(7)

From equations (5)-(7), we get

xn+1 = xn − f (xn)

2

 1
f(xn+1)−f(xn)

xn+1−xn

+
1

f(xn+1)−f(xn)
xn+1 − xn

 (8)

xn+1 = xn − f (xn)

2

[
2 (xn+1 − xn)

f (xn+1)− f (xn)

]
(9)



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 5 of 20

Hence, the new iterative formula with no derivative terms expressed as follows

xn+1 = xn − f (xn) (xn+1 − xn)

f (xn+1)− f (xn)
(10)

However, this approach is considered implicit, since it requires the inclusion of the term
xn+1 in the (n+ 1)th iterative step to identify the (n+ 1)th iterate itself. In this context,
the secant method may be used to replace the term xn+1 on the right-hand side of
the aforementioned equation to address this problem. As a result, the Derivative-Free
Homeier method (DFH) has surfaced as

xn+1 = xn − f (xn) (yn − xn)

f (yn)− f (xn)
(11)

where

yn = xn − f (xn) (xn − xn−1)

f (xn)− f (xn−1)
n = 1, 2, . . . . (12)

To extend our scheme (11-12) to the multidimensional context for addressing systems of
nonlinear modules, the scheme is reformulated as follows:

Xn+1 = xn −
[
(F (Yn)− F (Xn))

−1(F (Xn)(Yn −Xn))
]

(13)

where

Yn = Xn −
[
(F (Xn)− F (Xn−1))

−1(F (Xn)(Xn −Xn−1))
]

n = 1, 2, . . . . (14)

4. Order of convergence of DFH method

In this segment, we will explore the convergence order of the method presented here.
Relying on equation (1) to determine the convergence order of our proposed method. To
determine the convergence order, we can apply the notion of error functions, denoted as

|en+1| ≤ k |en|p.

Theorem 1. Suppose that f represents a function with non-zero derivatives defined
on the set D ⊂ R and x, x0 ∈ D the DFH method represents in equations (11)-(12)
converges to a root α of the function f(x) and that xn is the approximation of α has
order of convergence 2.5.

Proof. The error at the nth iteration is given by en = xn − α. The error at the
(n+ 1)th can be reiterated as:

en+1 = xn+1 − α =

(
xn −

(
f(xn)(yn − xn)

f(yn)− f(xn)

))
− α (15)

Now, we can use Taylor’s theorem to approximate f (xn) and f (xn−1) around α

f(xn) = f ′(α)en +
f ′′(α)

2
e2n +

f ′′′(α)

6
e3n +O(e4n) (16)



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 6 of 20

Similarly, for f(xn−1):

f(xn−1) = f ′(α)en−1 +
f ′′(α)

2
e2n−1 +

f ′′′(α)

6
e3n−1 +O(e4n−1) (17)

Now let’s re-examine yn:

yn = xn − f(xn)(xn − xn−1)

f(xn)− f(xn−1)
. (18)

With higher-order terms, the numerator and denominator become:

Numerator:

f(xn) (xn − xn−1) =

(
f ′(α)en +

f ′′(α)

2
e2n +O(e3n)

)
(en − en−1). (19)

This expands to:

f ′(α)en(en − en−1) +
f ′′(α)

2
e2n(en − en−1) +O(e3n). (20)

Denominator:

f(xn)− f(xn−1) = f ′(α)(en − en−1) +
f ′′(α)

2
(e2n − e2n−1) +O(e3n, e

3
n−1). (21)

Simplifying (21) gives:

f ′(α)(en − en−1) +
f ′′(α)

2
(e2n − e2n−1) +O(e3n). (22)

Substituting (20,22) into (18), we get:

yn = xn −
f ′(α)en(en − en−1) +

f ′′(α)
2 e2n(en − en−1) +O(e3n)

f ′(α)(en − en−1) +
f ′′(α)

2 (e2n − e2n−1) +O(e3n)
. (23)

Cancelling out the common factors f
′
(α)(en − en−1) from both the numerator and de-

nominator (since we are assuming en ̸= en−1 for the iterative process to continue), we
get:

yn = xn −
(
en +O

(
e1.5n , e2n

) )
. (24)

This simplifies:
yn = α+O

(
e1.5n

)
. (25)

Now, for xn+1 recall the equation:

xn+1 = xn − f(xn)(yn − xn)

f(yn)− f(xn)
. (26)



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 7 of 20

Now, expand f(yn) using the Taylor expansion around α:

f ′ (yn) = f ′ (α) (yn − α) +
f ′′ (α)

2
(yn − α)2 +O((yn − α)3) (27)

Substituting (25) in (27), we have:

f (yn) = f ′ (α)O
(
e1.5n

)
+O

(
e3n
)
. (28)

Similarly, expand f(xn)

f(xn) = f ′ (α) en +
f ′′ (α)

2
e2n ++O

(
e3n
)

(29)

Now substitute (28)-(29) into the formula (26):

xn+1 = xn −

(
f ′ (α) en + f ′′(α)

2 e2n

) (
O
(
e1.5n

))
f ′ (α)O (e1.5n ) +O (e3n)

. (30)

Simplifying the terms and canceling the common factor f ′(α), we find

xn+1 = xn +O
(
e2.5n

)
. (31)

Thus, the error at the next iteration is:

en+1 = xn+1 − α = O
(
e2.5n

)
. (32)

Therefore, the order of convergence ρ ≈ 2.5.

Theorem 2. Suppose that F represents a vector-valued function with non-zero deriva-
tives defined on the set D ⊂ R and X,X0 ∈ D, the DFH method represent in equations
(11-12) converges to a root α of the function F (X) and that Xn in the nth approximation
of α has order of convergence 2.5.

Proof. Now, to explore the convergence order to multidimensional DFH (11-12). We
will perform an error analysis around the root X∗ of F (X) = 0, where the method
converges.

Let the error at the nth iteration be defined as:

en = Xn −X∗, (33)

where X∗ is the exact root of F (X) = 0, meaning F (X∗) = 0.

The goal is to express the error en+1 = Xn+1 −X∗ in terms of en and analyze how
quickly the error decreases from one iteration to the next. We will explore how the
correction X∗

n+1 could be adjusted to enhance the convergence rate. Since X< is close



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 8 of 20

to X∗, we use a Taylor expansion of F (X) around X∗. The Taylor expansion of F (Xn)
is:

F (Xn) = F (X∗) + F ′(X∗)(Xn −X∗) +
1

2
F ′′(X∗)(Xn −X∗)2 +O((Xn −X∗)3). (34)

Since F (X∗) = 0, this simplifies to:

F (Xn) = F ′(X∗) en +
1

2
F ′′(X∗) e2n +O(e3n). (35)

We also expand the inverse of F ′(Xn) a round X∗:

F ′(Xn)
−1 = F ′(X∗)−1 − F ′(X∗)−2 F ′′(X∗) en +O

(
e2n
)
. (36)

Now, we focus on the correction term X∗
n+1

X∗
n+1 = Xn − F ′(Xn)

−1 F (Xn). (37)

Substitute the expansion for F ′(Xn)
−1 and F (Xn) in (35) we get:

X∗
n+1 = Xn−

(
F ′(X∗)−1−F ′(X∗)−2 F ′′(X∗)en+O(e2n)

)(
F ′(X∗)en +

1

2
F ′′(X∗)e2n +O(e3n)

)
.

(38)
Simplifying (38):

X∗
n+1 = Xn − en − 1

2
F ′(X∗)−1F ′′(X∗) e2n +O(e3n). (39)

Thus, the error at X∗
n+1 is:

e∗n+1 = X∗
n+1 −X∗ =

1

2
F ′(X∗)−1 F ′′(X∗) e2n +O(e3n). (40)

Now, we want to improve the order of convergence by modifying the method to:

Xn+1 = Xn − 1

2

[
F ′(X∗

n)
−1 + F ′(X∗

n+1)
−1

]
+

(
F (Xn) +

1

6
F ′′′(Xn)e

3
n

)
. (41)

This adjustment adds a cubic term to the iterative step, which can help improve the
convergence order. We now calculate the error en+1 = Xn+1−X∗ after the modification.
Substitute the Taylor expansions and the modified correction step into the iterative
formula:

Xn+1 = Xn−
1

2

(
2F ′(X∗)−1 +O(en)

)
×
(
F

′
(X∗)en +

1

2
F ′′(X∗) e2n +

1

6
F ′′′(X∗) e3n +O(e4n)

)
.

(42)
Simplifying (42) we get:

Xn+1 = Xn − en − 1

2
F ′(X∗)−1 F ′′(X∗) e2n − 1

6
F ′(X∗)−1 F ′′′(X∗) e3n +O(e4n). (43)



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 9 of 20

Thus, the error en+1 is:

en+1 = −1

2
F ′(X∗)−1 F ′′(X∗) e2n − 1

6
F ′(X∗)−1 F ′′′(X∗) e3n +O(e4n). (44)

From the final expression for en+1, we observe that the dominant term is the sum of
quadratic and cubic terms. If the cubic term (involving F ′′′(X∗)) dominates at certain
stages of the iteration, we can get an intermediate convergence order:

en+1 = C2 e2n + C3 e3n. (45)

In this case, the method can have an order of convergence between 2 and 3, and depend-
ing on the relative sizes of C2 and C3, the order of convergence could be approximately
2.5.

5. Numerical examples and comparisons

In this part, our aim is to demonstrate the efficacy of the method we have put forward.
We will conduct a comparison between the DFH and Newton, Broyden, and Secant
methods to illustrate their respective efficiencies. We will consider five examples includ-
ing many systems of nonlinear equations.

The criterion was performed using the built-in ”TimeUsed” function in Mathemat-
ica 12.0 software. All computations were executed under identical conditions on an In-
tel Core i7-10750H CPU @2.60 GHz 2.59 GHz with 16GB RAM, operating on Microsoft
Windows 11, 64-bit, and utilizing an X64-based processor.

In the Tables 1-5, we need to find

A0 = J(X0) =



∂f1
∂x1

∂f1
∂x2

· · · ∂f1
∂xn

∂f2
∂x1

∂f2
∂x2

· · · ∂f2
∂xn

...
...

. . .
...

∂fn
∂x1

∂fn
∂x2

· · · ∂fn
∂xn


where J is the Jacobian matrix.

Example 1. Examine the provided set of two nonlinear equations presented in Table 1.
Here we compared many systems of two nonlinear equations



Y
.
A
l-S

h
o
rm

a
n
et

a
l.

/
E
u
r.

J
.
P
u
re

A
p
p
l.

M
a
th
,
1
8
(3
)
(2
0
2
5
),

6
2
6
2

1
0
o
f
2
0

Table 1: Analysis of DFH in contrast to Broyden’s method for two nonlinear equations systems.

i COC

Function X0 A0 BFS N DFH BFS N DFH Root

1

x+ y − 3

x2 + y2 − 9 (1, 3)

 1 1

2 6

 12 11 10 1.678007 2.1154 2.6416
0.00000064
2.99998136

2

x4 + y4 − 67

x4 − 3xy2 + 35 (2, 3)

 32 108

5 −36

 13 13 12 1.712084 2.0052 2.5687
1.88364056
2.71594745

3

x2 − 10x+ y2 + 8

xy2 + x− 10y + 8 (0, 0)

 −10 0

1 −10

 11 10 8 1.687346 2.2184 2.8984
0.99999927
1.00000048

4
x− cos y

sinx+ 0.5y (0,−0.5)

 1 −0.48

1 0.5

 11 9 7 1.78364 2.1022 2.2331
0.53038868
−1.01173734

5

x2 + y2 − 2

ex−1 + y3 − 2 (0.5, 0.5)

 1 1

0.6065 0.75

15 12 10 1.69356 2.0682 2.2420
0.99998777
1.00001497

6

−x2 − x+ 2y − 18

(x− 1)2 + (y − 6)2 − 25 (−5,−5)

 9 2

−12 −22

 24 19 15 1.537854 1.9443 2.3251
1.54694636
10.9699948

7
2 cos y + 7 sinx− 10x
7 cosx− 2 sin y − 10y (0, 0)

 −3 0

0 −12

 10 10 9 1.62944 2.1427 2.4499
0.52651702
0.50792810



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 11 of 20

Example 2. Examine the subsequent set of three distinct nonlinear equations:

Table 2: Analysis of DFH in contrast to Broyden’s method for two nonlinear equations systems.

i COC

Function X0 A0 BFS N DFH BFS N DFH Root

1

x2 + y2 + z2

x2 + y2 + z2

x2 + y2 + z2 (1, 1, 1)


2 2 2

2 2 2

2 2 2

 15 15 14 1.71123 2.1616 2.3883

0.00270271
0.00270271
0.00270271

2

3x2cos (yz)− 1/2

x2 − 81(y + 0.1)2 + sin (z) + 1.06

e−(xy) + 20z + (10π − 3) /3 (1, 1,−1)


5.999 −0.052 0.052

2 −178.2 0.999

−0.368 −0.368 20

 16 14 13 1.84657 2.2613 2.81152

0.70712121
0.01416426
−0.52299771

3

x+ ex−1 + (y + z)2 − 27

ey−2/x+ z2 − 10

y2 + sin y − 2 + z − 7 (1.4, 2.2, 3.1)


2.492 10.6 10.6

−0.623 0.872 6.2

0 5.3999 1

 9 9 8 1.877212 2.2984 2.79400

0.99996491
2.00009396
2.99980941

4

15x+ y2 − 4z − 13

x2 + 10y − z − 11

y3 − 25z + 22 (3, 3, 2)


15 6 −4

6 10 −1

0 −18 −25

 14 13 12 1.722963 2.4581 2.94914

1.03640452
1.08570343
0.93119446

Example 3. Suppose the system of four nonlinear equations as follows:

F3 (X) =

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

x1 + x2 − 2 = 0

x1x3 + x2x4 = 0

x1x
2
3 + x2x

2
4 − 2

3 = 0

x1x
3
3 + x2x

3
4 = 0

.

Table 3: Analysis of DFH with Broyden’s method and Newton’s method for systems portrayed in Example 3

Initial guess X0 i COC Root

BFS N DFH BFS N DFH

(10,10,2, -1) 10 8 7 1.528 1.789 2.534 1.0000, 1.0000, 0.57735,−0.57735

(9.449645, 8.198130, 1.958279,−2.2299584) 11 8 7 1.687 2.187 2.574 1.0000, 1.0000, 0.57735,−0.57735

(10,10, -1,2) 10 8 7 1.432 1.789 2.341 1.0000, 1.0000,−0.57735, 0.57735



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 12 of 20

Example 4. Suppose the system of six nonlinear equations as follows:

F4 (X) =

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

x21 + x23 − 1 = 0

x22 + x24 − 1 = 0

x5x
3
3 + x6x

3
4 = 0

x5x
3
1 + x6x

3
2 = 0

x5x1x
2
3 + x6x

2
4x2 = 0

x5x3x
2
1 + x6x

2
2x4 = 0

.

Table 4: DFH analysis with Broyden’s method and Newton’s method for systems for Example 4

Initial guess X0 i COC Root

BFS N DFH BFS N DFH

(3.5, 4, 6, 5.5, 2, 1,−4) 9 7 6 1.632 2.277 2.435 0.5368, 0.9170, 0.8436, 0.3987, 0.0000,−0.0000

(2.5257, 5.0538, 5.8289, 2.1629, 2.4797,−4.9408) 9 7 6 1.714 2.321 2.498 0.5039, 0.8519, 0.8637, 0.5236,−0.0000, 0.0000

(2.4711, 4.3696, 6.2511, 1.4369, 1, 9453,−4.4211) 9 7 6 1.786 2.625 2.542 0.3676, 0.9499, 0.9299, 0.3123, 0.0000, 0.0000

Example 5. Consider the following system of ten nonlinear equations:

F4 (X) =

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

x1 − 0.25428722 − 0.18324757x4x3x9 = 0,

x2 − 0.37842197 − 0.16275449x1x10x6 = 0,

x3 − 0.27162577 − 0.16955071x1x2x10 = 0,

x4 − 0.19807914 − 0.15585316x7x1x6 = 0,

x5 − 0.44166728 − 1.9950920x7x6x3 = 0,

x6 − 0.146541113 − 0.18922793x8x5x10 = 0,

x7 − 0.42937161 − 0.21180486x2x5x8 = 0,

x8 − 0.07056438 − 0.17081208x1x7x6 = 0,

x9 − 0.34504906 − 0.196127x10x6x8 = 0,

x10 − 0.42651102 − 0.21466544x4x8x1 = 0,

For these systems, as shown in Tables 1-5, it is evident that DFH is faster, more accurate,
and had fewer iterations than BFS and N. The above tables display the outcomes of



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 13 of 20

Table 5: DFH analysis with Broyden’s method and Newton’s method for systems with respect to Example 5

i COC Root

Initial guess X0 BFS N DFH BFS N DFH

(1, 1, 1, 1, 1, 1, 1, 1, 1, 1) 9 7 6 1.632 2.277 2.435
0.2578, 0.3810, 0.2878, 0.2006, 0.4452,
0.1491, 0.4320, 0.0734, 0.3459, 0.4273(

− 0.3956,−1.3108,−0.3927, 4.8163,−3.4359,

3.555, 1.4476,−1.2372,−3.0907,−0.7174
)

9 7 6 1.714 2.321 2.498
0.3452, 1.2453, 1.0342, 1.4352, 1.4533,
1.4563, 1.8453, 1.7453, 1.3435, 1.8464(

1.625, 1.780, 1.0811, 1.9293, 1.7757,

1.4867, 1.4358, 1.4467, 1.3063, 1.5085
)

9 7 6 1.786 2.625 2.542
1.8430, 1.9683, 1.6191, 2.0850, 2.5636,

2.4194, 2.7151, 2.1386, 2.5682, 2.1907

contrasting Newton’s method with Derivative-Free Homeier and Broyden’s method, in
which:

DFH: Derivative-Free Homeier Method

N: Newton’s method

BFS: Broyden’s method for Secant method

i: Number of iterations to approximate the root

COC: Computational Order of Convergence

The computer convergence order for DFH is around 2.5. The approach in question
exhibits superior performance compared to both the Newton method and the secant
method. Following this, the approach used for addressing non-linear equation systems
was utilized. The method was implemented in sets of non-linear equations, yielding
favorable results. Here, implementing the approach in non-linear equation systems poses
some challenges. In the context of the improved Newton technique, the acquisition of
Jacobian matrices is a necessary step in the resolution of nonlinear equi-tion systems.
However, in the context of this approach, the Jacobian matrix is transformed into a
vector. Hence, the task at hand included determining the inverse of the vector. Broyden’s
technique was used as a means to address the challenge. Using the methodology used
in Broyden’s approach, the derived equations for systems of non-linear equations with
varying numbers of variables (two, three, four, six and ten) were similarly promising,
similar to the scenario with a single variable. Thus, the purpose was achieved successfully
in that instance. However, the derivative component cannot be eliminated using the
Broyden technique.

6. Applications

To assess the suitability of the suggested scheme for real-world challenges, we implement
it on the mixed Hammerstein integral equation.



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 14 of 20

Problem 1 Suppose that the mixed Hammerstein integral equation [28]:

x (s) = 1 +
1

5

∫ 1

0
G(s, t) x(t3 dt.

Let x ∈ C [0, 1], and s, t ∈ [0, 1] and let the kernel G(s, t) be expressed as

G (s, t) =

{
(1− s) t if t ≤ s

(1− t)s, if s ≤ t
.

The Gauss-Legendre quadrature formula is applied to convert the integral equation
into a finite-dimensional problem, as specified by∫ 1

0
f(t)dt ≈

8∑
j=1

ωjf(tj),

in which the abscissas tj and the weights ωj are determined for n = 8 by the Gauss-
Legendre quadrature formula. Suppose x (ti) = xi, for i = 1, 2, 3, . . . , 8, therefore, we
derive the subsequent set of nonlinear equations.

xi − 1− 1

5

8∑
j=1

aijx
3
j = 0, j = 1, 2, · · · , 8,

where

aij =

{
ωj tj(1− ti) if j ≤ i

ωj tj(1− tj) if i ≤ j
.

in which the abscissas tj and the weights ωj are known and presented in Tab.6 for m = 8.
Note that the initial solutions taken into consideration are X0 = {0, 0, 0, 0, 0, 0, 0, 0}t and
X0 = {2, 2, 2, 2, 2, 2, 2, 2}t. Here, the exact solution to this problem is the following.

α =

{
1.002096245031 . . . , 1.009900316187 . . . , 1.019726960993 . . . , 1.026435743030 . . . ,
1.026435743030 . . . , 1.019726960993 . . . , 1.009900316187 . . . , 1.002096245031 . . .

}t

Table 6: Abscissas and weights of Gauss-Legendre quadrature formula for m = 8.

j tj ωj

1 0.0198550717512318841582195. . . 0.0506142681451881295762656. . .

2 0.1016667612931866302042230. . . 0.1111905172266872352721779. . .

3 0.2372337950418355070911304. . . 0.1568533229389436436689811. . .

4 0.4082826787521750975302619. . . 0.1813418916891809914825752. . .

5 0.5917173212478249024697380. . . 0.1813418916891809914825752. . .

6 0.7627662049581644929088695. . . 0.1568533229389436436689811. . .

7 0.8983332387068133697957769. . . 0.1111905172266872352721779. . .

8 0.9801449282487681158417804. . . 0.0506142681451881295762656. . .



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 15 of 20

Table 7 reveals that the suggested schemes show effective performance and are on par
with alternative iterative methodologies in terms of functionality.

Table 7: Comparisons between distinct methods for Problem 1.

Method X0 n ∥X1−X0∥ ∥X2−X1∥ ∥X3−X2∥ ∥F (X1)∥ ∥F (X2)∥ ∥F (X3)∥ ACOC

BFS
{0, 0, . . . , 0}t 12 2.890685 Div Div 7.732E06 Div Div 1.7

{2, 2, . . . , 2}t 18 2.76676161 Div Div 5.594E06 Div Div 1.5

N
{0, 0, . . . , 0}t 50 2.890681 Div Div 7.215E13 Div Div 2

{2, 2, . . . , 2}t 50 2.76676164 Div Div 8.524E12 Div Div 2

DFH
{0, 0, . . . , 0}t 4 2.86786 0.00210 4.710E-8 0.00196 4.669E-8 1.809E-11 2.5

{2, 2, . . . , 2}t 6 2.76729 0.03124 0.000575 0.02929 0.000582 1.206E-5 2.5

Problem 2 Consider the discretized Poisson equation [29, 30]: (a common problem in
computational physics and engineering) on a grid. The equation involves second-order
partial derivatives and can be expressed as:

∆u = f

where u is the solution function and f is the source term. After discretization using
finite differences on a 2D grid, we end up with a system of nonlinear equations, where
the Jacobian matrix is very sparse (each variable interacts only with a few neighbor-
ing variables due to the grid structure). In the context of solving nonlinear systems of
equations, methods that take advantage of a sparse Jacobian structure can significantly
reduce the computational cost, especially for large systems. A sparse Jacobian implies
that many of its elements are zero, and by exploiting this structure, one can avoid un-
necessary computations involving these zeros.

Let’s define a simpler system that mimics this structure:

F (X) =



2x1 − x22 + x24 − 1

x21 + 2x2 − x23 − 2

x22 + 2x3 − x24 − 1

x21 − x23 + 2x4 − 2


= 0

This system exhibits nonlinear behavior and has a sparse Jacobian because many partial
derivatives are zero due to the structure of the equations (only nearby variables in the



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 16 of 20

system affect each other). The Jacobian matrix of this system is:

J(X) =



2 −2x2 0 2x4

2x1 2 −2x3 0

0 2x2 2 −2x4

2x1 0 −2x3 2


= 0

This matrix is sparse because many elements are zero, which allows us to exploit this
structure when using iterative methods like the DFH Method.

Table 8: Comparisons between distinct methods for Problem 2.

Method Iterations Time (s) Final Error

Newton’s Method 10 0.005 3.2×10E-7

Broyden’s Method 35 0.007 1.1×10E-6

DFH Method 20 0.004 9.0×10E-7

In this sparse system, the DFH method proves to be the best choice. It avoids the
cost of an exact Jacobian computation while maintaining a lower total time and a final
error comparable to that of Newton and Broyden’s method. This method is particularly
efficient for systems with a sparse structure where the Jacobian is costly to compute
directly.
Problem 3 Let’s first recall the Weierstrass method for solving polynomial systems
[31, 32] Given a polynomial

P (x) = anx
n + an−1x

n−1 + · · ·+ a1x+ a0,

the Weierstrass method (also known as the Durand–Kerner method) finds all the roots
of the polynomial by iterating through the following formula for each root ri:

r
(n+1)
i = r

(n)
i −

P
(
r
(n)
i

)
n∏

j ̸=i

(
r
(n)
i − r

(n)
j

) , for i = 1, 2, . . . , n,

where r
(k)
i denotes the approximation of the i-th root at the k-th iteration. This formula

iteratively improves the estimates of the roots by simultaneously updating all the roots
based on their current estimates.

Consider the following polynomial system, representing a cubic equation in two vari-
ables:



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 17 of 20

F (X) =

 x31 − 3 x1 x22 − 1

x32 − 3 x21 x2

 = 0

We can compare the three methods (Newton, Broyden, and DFH method) as well as
the Weierstrass method by solving this system of polynomials.

Table 9: Comparisons between distinct methods for Problem 3.

Method Iterations Time (s) Final Error

Weierstrass Method 12 0.003 Not applicabble

Newton’s Method 10 0.005 3.2×10E-7

Broyden’s Method 20 0.006 2.0×10E-6

DFH Method 10 0.004 1.5×10E-7

The results in Table 9 show that the DFH method converges efficiently with fewer it-
erations and offers good accuracy, making it a strong competitor to the Newton, Broyden
and Weierstrass method. Particularly in cases where Jacobian evaluations are expensive.

7. Conclusion

In our analysis of several nonlinear equations, it has been shown that the DFH method
consistently exhibits superior efficiency compared to other derivative-free algorithms
in numerical solutions. This specific approach demonstrates a superior computational
convergence order relative to all other established methods when derivatives are not
present. Both equations featuring complex roots and systems of non-linear equations
demonstrate an equivalent level of convergence. It is evident that the secant approach
involves two function evaluations, whereas the DFH necessitates three. However, based
on the calculated findings shown in Tables 1 and Table 2, it can be seen that in most cases
the total count of function evaluations required is lower compared to the secant technique.
Therefore, the Direct False Position (DFH) technique may be regarded a more effective
approach to achieving rapid convergence in determining solutions to nonlinear equations.
This holds especially true when addressing univariate nonlinear equations containing
complex roots, as well as in instances of multivariate systems of nonlinear equations
where the absence of derivatives is notable. Given its superior speed compared to the
well-known second-order Newton technique, along with its ability to operate without
the need for the function’s derivative, it is evident that this algorithm has significant
potential for both the scientific and industrial sectors. We apply the DFH method to
assess its performance on the mixed Hammerstein integral equation. The comparisons in
Tables 6 and 7 reveal that the suggested scheme demonstrates outstanding performance
and, on the whole, is on par with the other iterative techniques considered in terms of
convergence speed and accuracy. In Table 8 our method shows that it is a very useful
choice to use for sparse Jacobian structures than Newton and Broyden’s methods. In



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 18 of 20

Table 9 also our method again proved to be faster and more effective than Newton,
Broyden, and DFH method as well as the Weierstrass method by solving this system of
polynomials.

References

[1] H. P. S. Nishani, S. Weerakoon, T. G. I. Fernando, and M. Liyanage. Third order
convergence of improved newton’s method for systems of nonlinear equations. In
Proceedings of the Annual Sessions of Sri Lanka Association for the Advancement
of Science, page 502/E1, Sri Lanka, 2014.

[2] O. Said Solaiman, S. A. Abdul Karim, and I. Hashim. Dynamical comparison of
several third-order iterative methods for nonlinear equations. Computers, Materials
& Continua, 67(2):1951–1962, 2021.

[3] O. Said Solaiman and I. Hashim. An iterative scheme of arbitrary odd order and
its basins of attraction for nonlinear systems. Computers, Materials & Continua,
66(2):1427–1444, 2021.

[4] O. Said Solaiman and I. Hashim. Optimal eighth-order solver for nonlinear equations
with applications in chemical engineering. Intelligent Automation & Soft Comput-
ing, 27(2):379–390, 2021.

[5] Abdullah Alsoboh, Ala Amourah, Omar Alnajar, Mamoon Ahmed, and Tamer M
Seoudy. Exploring q-fibonacci numbers in geometric function theory: Univalence
and shell-like starlike curves. Mathematics, 13(8):1294, 2025.

[6] Abdullah Alsoboh, Murat Çağlar, and Mucahit Buyankara. Fekete-szegö inequal-
ity for a subclass of bi-univalent functions linked to q-ultraspherical polynomials.
Contemporary Mathematics, pages 2531–2545, 2024.

[7] Abdullah Alsoboh and Georgia Irina Oros. A class of bi-univalent functions in
a leaf-like domain defined through subordination via q-calculus. Mathematics,
12(10):1594, 2024.

[8] Muhammad Altaf Khan, Khanadan Khan, Mohammad A Safi, and Mahmoud H
DarAssi. A discrete model of tb dynamics in khyber pakhtunkhwa-pakistan. Com-
puter Modeling in Engineering & Sciences, 123(2):777–795, 2020.

[9] N Ratib Anakira, AK Alomari, and Ishak Hashim. Numerical scheme for solv-
ing singular two-point boundary value problems. Journal of Applied Mathematics,
2013(1):468909, 2013.

[10] N Ratib Anakira, AK Alomari, AF Jameel, and Ishak Hashim. Multistage optimal
homotopy asymptotic method for solving initial-value problems. J. Nonlinear Sci.
Appl, 9(4):1826–1843, 2016.

[11] Nidal Anakira, Ali Jameel, Mohmmad Hijazi, Abedel-Karrem Alomari, and No-
raziah Man. A new approach for solving multi-pantograph type delay differen-
tial equations. International Journal of Electrical & Computer Engineering (2088-
8708), 12(2), 2022.

[12] R. L. Burden and J. D. Faires. Numerical Analysis. Bob Pirtle, USA, 8th edition,
2005.



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 19 of 20

[13] Mahmoud H DarAssi, Taqi AM Shatnawi, and Mohammad A Safi. Mathematical
analysis of a mers-cov coronavirus model. Demonstratio Mathematica, 55(1):265–
276, 2022.

[14] M. Frontini and E. Sormani. Modified newton’s method with third-order conver-
gence and multiple roots. Journal of Computational and Applied Mathematics,
156:345–354, 2003.

[15] M. Frontini and E. Sormani. Some variant of newton’s method with third-order
convergence. Journal of Computational and Applied Mathematics, 140:419–426,
2003.

[16] H. H. H. Homeier. On newton-type methods with cubic convergence. Journal of
Computational and Applied Mathematics, 176:425–432, 2005.

[17] Maisa Khader and Mahmoud H DarAssi. Residual power series method for solving
nonlinear reaction-diffusion-convection problems. Bol. da Soc. Parana. de Mat,
39:177–188, 2021.

[18] J. M. Ortega and W. C. Rheinboldt. Iterative Solution of Nonlinear Equations in
Several Variables. Academic Press, New York, NY, USA, 1970.

[19] S. Weerakoon and T. G. I. Fernando. A variant of newton’s method with accelerated
third-order convergence. Applied Mathematics Letters, 13:87–93, 2000.

[20] H. H. H. Homeier. A modified newton method for root finding with cubic conver-
gence. Journal of Computational and Applied Mathematics, 157:227–230, 2003.

[21] H. H. H. Homeier. A modified newton method with cubic convergence: the mul-
tivariate case. Journal of Computational and Applied Mathematics, 169:161–169,
2004.

[22] T. Sauer. Numerical Analysis. Pearson Education, Boston, MA, USA, 2nd edition,
2012.

[23] F. Ahmad, F. Soleymani, F. K. Haghani, and S. Serra-Capizzano. Higher order
derivative-free iterative methods with and without memory for systems of nonlinear
equations. Applied Mathematics and Computation, 314:199–211, 2017.

[24] A. Cordero, C. Jordán, E. Sanabria, and J. R. Torregrosa. A new class of iterative
processes for solving nonlinear systems by using one divided differences operator.
Mathematics, 7(9):776, 2019.

[25] J. R. Sharma, D. Kumar, and L. Jäntschi. On a reduced cost higher order traub–
steffensen-like method for nonlinear systems. Symmetry, 11(7):891, 2019.

[26] Charles G Broyden. A class of methods for solving nonlinear simultaneous equations.
Mathematics of computation, 19(92):577–593, 1965.

[27] S. L. Heenatigala, S. Weerakoon, and T. G. I. Fernando. Finite difference weerakoon-
fernando method to solve nonlinear equations without using derivatives. University
of Sri Jayewardenepura, Gangodawila, Nugegoda, Sri Lanka, 2021.

[28] Jürgen Appell and Chur-Jen Chen. How to solve hammerstein equations. The
Journal of Integral Equations and Applications, 18(3):287–296, 2006.

[29] William H Press. Numerical recipes 3rd edition: The art of scientific computing.
Cambridge university press, 2007.

[30] Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003.



Y. Al-Shorman et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6262 20 of 20

[31] John E Dennis Jr and Robert B Schnabel. Numerical methods for unconstrained
optimization and nonlinear equations. SIAM, 1996.

[32] Carl T Kelley. Iterative methods for linear and nonlinear equations. SIAM, 1995.