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

Sixteenth-Order Steffensen-Ostrowski Approach for Nonlinear
Problems with Applications in Celestial, Predator-Prey and

Neural Activation
Saima Akram1,2, Muhammad Bilal Riaz3,5,*, Hareem Khalid2, Fouzia Amir*4,
Faiza Akram2, Dexkanov Suxrob Sobirovich4

1 Department of Mathematics, Government College Women University Faisalabad, Faisalabad 38000
Pakistan
2 Center for Advanced Studies in Pure and Applied Mathematics, Bahauddin Zakariya University,
Multan, 60000, Pakistan
3 IT4Innovations, VSB – Technical University of Ostrava, Ostrava, Czech Republic
4 Centre for Research and Innovation, Asia International University, Bukhara, 200100, Uzbekistan
5 Applied Science Research Center, Applied Science Private University, Amman, Jordan

Abstract. The increasing demand for accurate and efficient solutions to nonlinear equations, driven by ad-
vancements across diverse research and engineering fields, highlights the critical need for innovative com-
putational methods. This study addresses this need by introducing a novel derivative-free sixteenth-order
iterative scheme derived from a weighted Steffensen-Ostrowski-type family. According to the Kung-Traub
conjecture, this scheme is designed to achieve optimal convergence using only five function evaluations per
iteration. A key innovation lies in employing a bivariate weight function in the third step and Lagrange
interpolation in the fourth step, ensuring high accuracy and computational efficiency by avoiding derivative
evaluation. The extensive convergence analysis shows that the proposed scheme is of sixteenth order and
is validated through applications to real-world problems, including Kepler’s celestial motion, an ideally
mixed reactor, predator-prey models, neural activation dynamics, and periodic ecosystem growth. Numer-
ical results demonstrate the superiority of the proposed scheme over existing four-point iterative schemes,
particularly in terms of absolute error and computational convergence order. Furthermore, graphical analy-
sis of complex polynomials illustrates the algorithm’s attraction basins, offering a wide range of choosing
from initial guesses to converge to the specified root more efficiently without divergence and hence showed
more stable behavior. In this way we achieved significant computational improvements and higher con-
vergence order. As the proposed scheme fall under the category of derivative free schemes, so it is more
general as compared to existing schemes in literature and is considered to be a good alternative to the
existing schemes especially where derivatives are unavailable.

2020 Mathematics Subject Classifications: 65H10, 65B99, 92B20

∗Corresponding author.
∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5822

Email addresses: saimaakram@gcwuf.edu.pk (S. Akram),
muhammad.bilal.riaz@vsb.cz (M. Riaz), faiza4akram@gmail.com (F. Akram),
fouziaur.abbasi@oxu.uz (F. Amir), hareemkhalid65@gmail.com (M. Khalid),
s.dekhanov@oxu.uz (S. Dekhanov)

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



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 2 of 26

Key Words and Phrases: Nonlinear equations, Iterative methods, Lagrange interpolation, Optimal conver-
gence order, Efficiency index, Basins of attraction

1. Introduction

The branch of computational mathematics that provides approximated solutions to a variety of
mathematical problems is known as Numerical Analysis. These problems arise from practical ap-
plications in algebra, geometry, calculus, and span fields such as natural sciences, social sciences,
engineering, medicine, and business. They often involve variables that change continuously. One
of the key challenges in these fields is solving nonlinear equations, defined as:

g(k) = 0, (1)

plays a significant role in computational and applied mathematics. Here, g : D ⊆ C → C is a
sufficiently differentiable function on the interval D. Nonlinear equations in mathematics come
in various forms, such as transcendental, integral, algebraic, and both ordinary and partial differ-
ential equations, often in combination. Analytic methods cannot solve these types of equations.
Approximating the roots is a standard technique in numerical analysis. Numerical methods are key
to tackling the complexity of nonlinear equations. A wealth of literature exists on solving these
equations, including foundational works by Ostrowski [1], Kung and Traub [2], and Petković [3].

There are two widely used techniques for estimating the roots of a nonlinear algebraic equa-
tion: one-point and multipoint iteration. Newton’s method was first presented by Thomas Simpson
in 1740 as an iterative method of solving nonlinear equations. Newton’s method, known for its
one-point iterative approach, exhibits quadratic convergence. However, it has limitations, particu-
larly its dependence on derivative evaluation and the necessity for a suitable initial guess to achieve
root convergence. A recent suggestion to avoid introducing extra functions such as evaluating first
or second derivatives, is to eliminate derivatives from the iteration process. As an example, this ap-
proach leads to the familiar Steffensen’s method [2], where the forward difference approximation
g(k j+g(k j))−g(k j)

g(k j)
substitutes the first-order derivative g′(kn) in Newton’s method as follow:

k j+1 = k j −
g(k j)

g[k j,v j]
,

where v j = k j + g(k j) and g[k j,v j] =
g(v j)−g(k j)

v j−k j
is the divided difference of first order. Both

the Steffensen approach and the Newton method acquire two function evaluations and converge
quadratically. Kung and Traub’s conjecture [2], posits that any without memory multipoint iter-
ative method has convergence order cannot surpass the upper limit of 2n−1, limited to n function
evaluations. With advances in digital computers, arithmetic, and symbolic computing, higher-
order multipoint methods have gained popularity. These methods offer more accurate and efficient
root estimates in fewer iterations, with an efficiency index [4] superior to that of Newton’s method.
Techniques such as weight functions, Taylor expansions, and dynamical analysis have played a key
role in developing both optimal and non-optimal multipoint methods with convergence orders up
to eight [5–17]. Despite these advancements, achieving optimal four-point methods remains a



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 3 of 26

significant challenge and an area of ongoing research. Some optimal sixteenth-order methods for
solving nonlinear equations have been developed by several researchers (for example see, [18–
25]). While most of these methods achieve high convergence orders, they often rely on derivative
evaluations, which are computationally intensive and not always feasible for all problems.

These challanges inspired us to create a novel four-point optimal iterative approach based
on the weighted family of Steffensen-Ostrowski type. It requires five functional evaluations and
gives sixteenth-order convergence per iteration. The first two step of the suggested method is
Ostrowski’s method and third step is constructed by using the bivariate weight function. For
the construction of the fourth step, Lagrange interpolation technique has been used. Lagrange
interpolation is highly advantageous for its simplicity, derivative-free formulation, and flexibility
in node selection, making it efficient for improving iterative methods for nonlinear equations.

This work aims to satisfy the Kung and Traub conjecture by achieving optimal sixteenth order
convergence supported by an extensive convergence analysis. To evaluate its accuracy and effec-
tiveness, the proposed approach is tested on real-world problems, including Kepler’s equation of
celestial motion, an ideally mixed reactor, predator-prey models, neural activation dynamics, and
periodic ecosystem growth. Numerical results demonstrate the superiority of the proposed scheme
over existing four-point iterative methods, particularly in terms of absolute error and convergence
order. Additionally, graphical analysis of complex polynomials shows the algorithm’s attraction
basins. This analysis illustrating that it converges efficiently to the specified root from a wide range
of initial guesses, with improved stability and minimum divergence. The proposed derivative-free
scheme is more general than existing methods in the literature and offers a promising alternative
in the cases where derivatives does not exist.

The rest of the paper is designed as: The formulation of a four-point iterative method along
with the analysis of convergence is presented in Section 2. Some particular cases of weight func-
tions are discussed in Section 3. Standard test functions and numerical experimentation along with
the comparison of the developed methods with existing methods of equivalent order are depicted
in Section 4. In Section 5, a detailed dynamical analysis of the presented approaches in a complex
plane is demonstrated using a graphical tool basins of attraction. Finally, concluding remarks are
given in Section 6.

2. A Formulation of Optimal sixteenth-order Scheme Using Lagrange Interpolation

In this section, an optimal iterative scheme has been presented which is based on Steffensen-
Ostrowski’s type weighted family.

2.1. Formulation of Scheme

For the construction of the scheme, take into account the three-point optimal eighth-order
iterative method proposed by Kanwar et al. [26], based on Ostrowski’s method:

y j = k j −
g(k j)

g[ζ j,k j]
, ζ j = k j +βg(k j)

3,

z j = y j −
g(y j)

2g[y j,k j]−g[ζ j,k j]
,



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 4 of 26

k j+1 = z j −
g(z j)

g[y j,z j]+g[ζ j,y j,z j](z j − y j)
E(λ ,µ), (2)

where, λ =
g(z j)
g(y j)

and µ =
g(y j)
g(k j)

. β is a parameter and β ∈ R\{0}.
With the help of Newton’s technique at the fourth step of the scheme (2) we obtain:

y j = k j −
g(k j)

g[ζ j,k j]
, ζ j = k j +βg(k j)

4,

z j = y j −
g(y j)

2g[y j,k j]−g[ζ j,k j]
,

t j = z j −
g(z j)

g[y j,z j]+g[ζ j,y j,z j](z j − y j)
E(λ ,µ),

k j+1 = t j −
g(t j)

g′(t j)
, (3)

where, β ∈ R\{0}, λ =
g(z j)
g(y j)

and µ =
g(y j)
g(k j)

.

The above scheme (3) does not meet the criteria for optimality as per the Kung-Traub conjec-
ture [2], since per iteration it necessitates six function evaluations. Consequently, introducing a
suitable approximation, we sought to reduce the number of function evaluations. The below given
approximation is employed to replace the derivative g′(t j):

g′(t j) =
(t j − z j)(t j − k j)(t j − y j)g(ζ j)

(ζ j − t j)(ζ j − k j)(ζ j − z j)(ζ j − y j)

+
(t j − z j)(t j −ζ j)(t j − k j)g(y j)

(y j − t j)(y j −ζ j)(y j − z j)(y j − k j)

+
(t j − z j)(t j −ζ j)(t j − y j)g(k j)

(k j − t j)(k j −ζ j)(k j − y j)(k j − z j)

+
(t j −ζ j)(t j − k j)(t j − y j)g(z j)

(z j − t j)(z j −ζ j)(z j − k j)(z j − y j)

+
g(t j)

(t j −ζ j)
+

g(t j)

(t j − k j)
+

g(t j)

(t j − y j)
+

g(t j)

(t j − z j)
. (4)

where, ζ j = k j +βg(k j)
4 and β ∈ R \ {0}. g′(t j) is the fourth-degree Interpolation by Lagrange

[27] that interpolates t j, z j, y j, k j, and ζ j and β ∈ R\{0}.
We formulate a four-step iterative method of order sixteenth utilizing a bivariate weight func-

tion and at each iterative step employing five evaluations of function. This formula achieved by
integrating the provided approximation (4) in the fourth step of the iterative scheme (3), as shown
below:

y j = k j −
g(k j)

g[ζ j,k j]
, ζ j = k j +βg(k j)

4,

z j = y j −
g(y j)

2g[y j,k j]−g[ζ j,k j]
,



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 5 of 26

t j = z j −
g(z j)

g[y j,z j]+g[ζ j,y j,z j](z j − y j)
E(λ ,µ),

k j+1 = t j −
g(t j)

ML
. (5)

where, β ∈ R\{0}, λ =
g(z j)
g(y j)

, µ =
g(y j)
g(k j)

, and

ML =
(t j − z j)(t j − k j)(t j − y j)g(ζ j)

(ζ j − t j)(ζ j − k j)(ζ j − z j)(ζ j − y j)

+
(t j − z j)(t j −ζ j)(t j − k j)g(y j)

(y j − z j)(y j − t j)(y j −ζ j)(y j − k j)

+
(t j − z j)(t j −ζ j)(t j − y j)g(k j)

(k j − y j)(k j − t j)(k j −ζ j)(ζ j − z j)

+
(t j −ζ j)(t j − k j)(t j − y j)g(z j)

(z j − t j)(z j −ζ j)(z j − k j)(z j − y j)

+
g(t j)

(t j −ζ j)
+

g(t j)

(t j − k j)
+

g(t j)

(t j − y j)
+

g(t j)

(t j − z j)
.

Consider an analytic function E : C2 → C in the vicinity of (0,0), with λ =
g(z j)
g(y j)

= O(e2
j)

and µ =
g(y j)
g(k j)

= O(e j). The following theorem demonstrates that the order of convergence of the
aforementioned scheme is optimal sixteenth-order.

Theorem 1. Let g be a sufficiently differentiable function defined on an open interval D, with
γ ∈ D being a simple zero of g. Assume that the initial guess k∗ is sufficiently close to γ ∈ D. Then,
the four-point iterative scheme described in (5) achieves an optimal convergence order of sixteen,
provided the following conditions for the weight function are satisfied:

E00 = 1 = E11,

E01 = 0 = E10 = E02,

E03 = −6. (6)

where Ei, j = [ ∂E(λ ,ρ)
∂λ iρ j

1
i! j! ](0,0) and i, j = 0, 1, 2, 3. It possess the error relation as follows:

k j+1 − γ = (
1
6

c3
3c9

2g−6c12
2 c4 − . . .− 1

6
gc3c13

2 b2 +3c3
3c9

2a)e16
j +O(e17

j ). (7)

where e j = k j − γ and c j =
g( j)(γ)
j!g′(γ) , j = 1, 2, 3 · · · .

Proof. Expanding g(k j) around the simple zero γ, using Taylor’s expansion, considering that
g(γ) = 0:

g(k j) = g′(γ)(e j + c2e2
j + c3e3

j + c4e4
j + c5e5

j + · · ·+ c16e16
j )+O(e17

j ). (8)



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 6 of 26

where c j =
g( j)(γ)
j!g′(γ) , j = 1, 2, 3 · · · . We define ζ j as :

ζ j = k j +βg(k j)
4,

ζ j − γ = e j +βe4
j +4βc2e5

j +2β (3c2
2 +2c3)e6

j + · · · . (9)

Now we expand g(ζ j) using Taylor expansion about γ as:

g(ζ j) = g′(γ)(e j + c2e2
j + c3e3

j +(c4 +β )e4
j +(c5 +6βc2)e5

j + · · ·). (10)

To compute divided difference g[ζ j,k j], we utilized equations (8) and (10) as follows:

g[ζ j,k j] =
g(ζ j)−g(k j)

w j − k j
,

= 1+2c2e j +3c3e2
j +4c4e3

j +(5c5 +βc2)e4
j + · · · . (11)

Now, equations (8) and (11) substituting in the first step of method (5) and we get:

y j − γ = c2e2
j +(2c3 −2c2

2)e
3
j +(3c4 −7c2c3 +4c3

2)e
4
j + · · · . (12)

Hence, equation (12) is achieved second order convergence of method (5) which is optimal. Taylor
expansion of g(yn) about γ as well as considering equation (12):

g(y j) = g′(γ)(c2e2
j +(2c3 −2c2

2)e
3
j +(3c4 −7c2c3 +5c3

2)e
4
j + · · ·). (13)

To determine the divided difference g[y j,k j], we used equations ( 8) and (13 ) subsequently, we
arrive at the expression as follows:

g[y j,k j] =
g(k j)−g(y j)

k j − y j
, (14)

= 1+ c2e j +(c1 + c2
2 + c3)e2

j +(c4 +3c3c2 −2c3
2)e

3
j + · · · .

In the second step of iterative method (5), we used equations (11), (13), and (14) to achieve fourth
order convergence:

z j − γ = (−c2c3 + c3
2)e

4
j +(−2c2c4 −4c4

2 +8c2
2c3 −2c2

3)e
5
j + · · · . (15)

Expanding g(z j) about γ by using of equation (15), we have:

g(z j) = g′(γ)(c2(−c3 + c2
2)e

4
j +(−2c2c4 −4c4

2 +8c2
2c3 −2c2

3)e
5
j + · · ·). (16)

Next, we derive the Taylor expansion of the divided differences used in the third step of the scheme
(5). Employing equations (13) and (16) we find g[y j,z j], as follows:

g[y j,z j] =
g(z j)−g(y j)

z j − y j
,

= 1+ c2
2e2

j −2c2(−c3 + c2
2)e

3
j + · · · . (17)



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 7 of 26

To find g[ζ j,y j], we used equations (10) and (13) and get the following expression:

g[ζ j,y j] =
g(y j)−g(ζ j)

y j −ζ j
,

= 1+ c2e j +(c3 + c2
2)e

2
j +(c4 −2c3

2 +3c2c3)e3
j + · · · . (18)

we calculated divided difference g[ζ j,y j,z j], by using equations (17) and (18):

g[ζ j,y j,z j] =
g[y j,z j]−g[ζ j,y j]

z j −ζ j
,

= c2 + c3e j +(c4 + c2c3)e2
j + · · · . (19)

Furthermore, we introduce two variables λ and µ , defined as follows:

λ =
g(z j)

g(y j)
= (c2

1 + c2
2 − c3)e2

j +(−2c4 +4c2c3 −2c3
2)e

3
j + · · · . (20)

and

µ =
g(y j)

g(k j)
= c2e j +(2c3 −3c2

2)e
2
j +(3c4 −10c2c3 +8c3

2)e
3
j + · · · . (21)

From equations (20) and (21) it was observed that λ exhibits second-order convergence, while µ

demonstrates first-order convergence. Now, the bivariate weight function E(λ ,µ) was expanded
using a Taylor series expansion up to the fifth term, as illustrated:

E(λ ,µ) = E00 +(λE10 +µE01)+
1
2!
(λ 2E20 +2λ µE11 +µ

2E02)

+
1
3!
(λ 3E30 +3λ

2
µE21 +3λ µ

2E12 +µ
3E03)+

1
4!
(λ 4E40

+4λ
3
µE31 +6λ

2
µ

2E22 +4λ µ
3E13 +λ

4E04)+
1
5!
(λ 5E50

+5λ
4
µE41 +10λ

3
µ

2E32 +10λ
2
µ

3E23 +5λ µ
4E14

+µ
5E05). (22)

In third step of scheme (5), we substituted equations (16), (17), (19), and (22) as:

t j − γ =
(
−c2c3 + c3

2 + c2E00c3 − c3
2E00

)
e4

j +(−4c4
2 +8c2

2c3 −2c2c4

−2c2
3 − c2

2(c2 − c3)E01(−4c4
2 +8c2

2c3 −2c2c4 −2c2
3)E00)e5

j

+
(
10c5

2 −30c3
2c3 −βc2

2 +12c2
2c4 +18c2c2

3 −3c2c5 −7c3c4
)

−1
2

E02c5
2 +

1
2

E02c3
2c3 −E10c5

2 +2E10c3
2c3 −E10c2c2

3 +7E01c5
2

−13E01c3
2c3 +4E01c2c2

3 −10E00c5
2 +30E00c3

2c3 +E00βc2
2

−12E00c2
2c4 −18E00c2c2

3 −3E00c2c5 +7E00c3c4

)
e6

j + · · · . (23)



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 8 of 26

Equation (23) demonstrates fourth-order convergence. To achieve eighth-order convergence, we
employed the following values:

E00 = 1 = E11,

E01 = 0 = E10 = E02,

E03 = −6. (24)

By substituting the values (24) in equation (23) we get the error term as:

t j − γ = (−c2
2c4c3 −

3
2
(b2c2

3c3
2)−

1
2
(ac2

3c3
2)+ · · ·− 1

2
(b2c7

2))e
8
j + · · · . (25)

By utilizing equation (25), we successfully attained optimal eighth-order convergence. Employing
Taylor series expansion, we computed g(t j), yielding the following expression:

g(t j) = g′(γ)(− 1
24

{c2(24c3c4c2 +12c6
2b2 −48c2

3c2
2 + c6

2g+12c6
2a−·· ·)}e8

j + · · ·). (26)

Firstly, we calculate (t j−k j)(t j−z j)(t j−y j)g(ζ j)
(ζ j−t j)(ζ j−k j)(ζ j−z j)(ζ j−y j)

term by utilizing equations (9), (10), (12), (15), and
(25) as:

(t j − k j)(t j − z j)(t j − y j)g(ζ j)

(ζ j − t j)(ζ j − k j)(ζ j − z j)(ζ j − y j)
=

−c2
2(−c3 + c2

2)

β
e j + · · · . (27)

Now, to determine the term (t j−ζ j)(t j−z j)(t j−y j)g(k j)
(k j−t j)(k j−y j)(k j−ζ j)(ζ j−z j)

we used (8), (9), (12), (15), and (25) given as:

(t j − k j)(t j − z j)(t j − y j)g(ζ j)

(ζ j − z j)(ζ j − k j)(ζ j − t j)(ζ j − y j)
=

−c2
2(−c3 + c2

2)

β
e j (28)

+
2c2(−7c2

2c3 +4c4
2 +2c2

3 + c2c4)

β
e2

j + · · · .

Also, to find (t j−z j)(t j−ζ j)(t j−k j)g(y j)
(y j−z j)(y j−t j)(y j−ζ j)(y j−k j)

we used (9), (12), (13), (15), and (25) as:

(t j − z j)(t j −ζ j)(t j − k j)g(y j)

(y j − z j)(y j − t j)(y j −ζ j)(y j − k j)
= (c3 − c2

2)e
2
j +(2c4 −2c2c3)e3

j + · · · . (29)

We used (9), (12), (15), (16), and (15) to find (t j−ζ j)(t j−k j)(t j−y j)g(z j)
(z j−t j)(z j−ζ j)(z j−k j)(z j−y j)

term given as:

(t j −ζ j)(t j − y j)(t j − k j)g(z j)

(z j − t j)(z j − k j)(z j − y j)(z j −ζ j)
= 1+(−c3 + c2

2)e
2
j +(−2c4 +2c2c3)e3

j + · · · . (30)

We computed g(t j)
(t j−w j)

by using equations (10), (25), and (26) as:

g(t j)

(t j −ζ j)
=

1
24

c2(24c2c3c4 −·· ·−36b2c3c4
2)e

7
j + · · · . (31)



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 9 of 26

Using equations (25) and (26) we evaluated g(t j)
(t j−k j)

as:

g(t j)

(t j − k j)
=

1
24

c2(24c2c3c4 −·· ·−36b2c3c4
2)e

7
j + · · · . (32)

g(t j)
(t j−y j)

is calculated by using equations (12), (25), and (26):

g(t j)

(t j − y j)
= (c2c3c4 −2c2

2c2
3 −·· ·− 3

2
b2c3c4

2)e
6
j + · · · . (33)

Last term g(t j)
(t j−z j)

is calculated by using equations (15), (25), and (26):

g(t j)

(t j − z j)
= (

1
24

c4
2g+

1
2

c4
2b2 + · · ·+ 1

2
b2c2

3)e
4
j + · · · . (34)

Substituting equations (27)-(34) in equation in ML, we obtained:

ML = 1− 1
12

c2
2(−12b2c3

3 +36b2c2
2c2

3 −·· ·−24ac3c4
2)e

8
j + · · · . (35)

Now, utilizing equations (26) and (35) in fourth step of iterative scheme (5) we obtained error term
as:

e j+1 = (
1
6
(c3

3gc9
2)−6c12

2 c4 − . . .− 1
6
(b2c13

2 gc3)+3c3
3ac9

2)e
16
j +O(e17

j ). (36)

Above error relation (36) demonstrates that the sixteenth-order convergence approach (5) is opti-
mal.

3. Some Particular Cases

We introduce two specific weight functions: Case 1 and Case 2.
Case 1: To obtain sixteenth-order convergence we utilize the bivariate polynomial of the fol-

lowing form:

H(λ ,µ) = 1+λ µ +5λ
2 −µ

3 −8µ
2
λ .

In the scheme’s third step (5), SFHM-1 is indicated and explained as follows:

y j = k j −
g(k j)

g[ζ j,k j]
, ζ j = k j +βg(k j)

4,

z j = y j −
g(y j)

2g[y j,k j]−g[ζ j,k j]
,

t j = z j −
g(z j)

g[y j,z j]+g[ζ j,y j,z j](z j − y j)
(1+λ µ +5λ

2 −µ
3 −8µ

2
λ +λ

3),

k j+1 = t j −
g(t j)

ML
. (37)



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 10 of 26

where β ∈ R\{0}, λ =
g(z j)
g(y j)

, µ =
g(y j)
g(k j)

and

ML =
(t j − z j)(t j − k j)(t j − y j)g(ζ j)

(ζ j − k j)(ζ j − t j)(ζ j − z j)(ζ j − y j)

+
(t j − z j)(t j −ζ j)(t j − y j)g(k j)

(k j − t j)(k j −ζ j)(k j − y j)(ζ j − z j)

+
(t j − z j)(t j −ζ j)(t j − k j)g(y j)

(y j − z j)(y j − t j)(y j −ζ j)(y j − k j)

+
(t j − k j)(t j −ζ j)(t j − y j)g(z j)

(z j − t j)(z j −ζ j)(z j − k j)(z j − y j)

+
g(t j)

(t j −ζ j)
+

g(t j)

(t j − k j)
+

g(t j)

(t j − y j)
+

g(t j)

(t j − z j)
.

Case 2: Using weight function of the form

H(λ ,µ) = 1+λ µ +(λ µ)3 −µ
3 −8µ

2
λ .

The new scheme, termed SFHM-2, takes on the following form in the third step of (5).

y j = k j −
g(k j)

g[ζ j,k j]
,ζ j = k j +βg(k j)

4,

z j = y j −
g(y j)

2g[y j,k j]−g[ζ j,k j]
,

t j = z j −
g(z j)

g[ζ j,y j,z j](z j − y j)+g[y j,z j]
(1+λ µ +(λ µ)3 −µ

3 −8µ
2
λ ),

k j+1 = t j −
g(t j)

ML
, (38)

where β ∈ R\{0}, λ =
g(z j)
g(y j)

, µ =
g(y j)
g(k j)

and

ML =
(t j − z j)(t j − k j)(t j − y j)g(ζ j)

(ζ j − t j)(ζ j − z j)(ζ j − k j)(ζ j − y j)

+
(t j − z j)(t j −ζ j)(t j − y j)g(k j)

(k j − y j)(k j − t j)(k j −ζ j)(ζ j − z j)

+
(t j − z j)(t j −ζ j)(t j − k j)g(y j)

(y j − z j)(y j − t j)(y j −ζ j)(y j − k j)

+
(t j − k j)(t j −ζ j)(t j − y j)g(z j)

(z j − t j)(z j −ζ j)(z j − k j)(z j − y j)

+
g(t j)

(t j −ζ j)
+

g(t j)

(t j − k j)
+

g(t j)

(t j − y j)
+

g(t j)

(t j − z j)
.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 11 of 26

4. Computational Analysis

In this section, we delve into examining the convergence behavior, effectiveness, robustness,
and validity of the schemes proposed in section 2. To assess accuracy, we employ the absolute
error between consecutive iterations |k j+1−k j| for the initial three iterations. The theoretical order
of convergence is confirmed through the calculation of the computational order of convergence
(COC).

p =
ln |(k j − k j−1)/(k j−1 − k j−2)|

ln |(k j−1 − k j−2)/(k j−2 − k j−3)|
.

The Maple 16 programming package uses multi-precision arithmetic with 2000 significant decimal
digits to ensure high accuracy and prevent the loss of significant digits. Different standard test
functions and some real-world problems have been used to check the efficiency of our proposed
schemes. The sixteenth-order convergence techniques that were employed for comparison are
listed below.

For the sake of comparison, we adopt the iterative scheme introduced by Sharma et al. [21],
referred to as CM-1. This method was constructed by using Rational Interpolation, outlined as:

z j = k j −
g(k j)

g′(k j)
,

w j = z j −
g′(k j)g(z j)

(g[z j,k j])2 ,

u j = w j −
1

h1 −h2 +h3

g′(k j)g(w j)

(g[w j,k j])2 ,

k j+1 = u j −
1

H1 −H2 +H3 −H4

g′(k j)g(u j)

(g[u j,k j])2 , (39)

where

h1 =
z j −w j

z j − k j
, h2 =

(w j − k j)
2g′(k j)

(z j −w j)(z j − k j)g[z j,k j]
,

h3 =
(w j − k j)g′(k j)

(z j −w j)g[w j,k j]
, H1 =

(z j −u j)(w j −u j)

(z j − k j)(w j − k j)
,

H2 =
(u j − k j)

2(u j −w j)g′(k j)

(z j − k j)(z j −u j)(z j −w j)g[z j,k j]
,

H3 =
(u j − k j)

2(z j −u j)g′(k j)

(w j −u j)(w j − z j)(w j − k j)g[w j,k j]
,

H4 =
(u j − k j)(2u j − z j −w j)g′(k j)

(w j −u j)(z j −u j)g[u j,k j]
.

We juxtapose the numerical results of our newly developed methods with an iterative four-
point scheme devised by Sharifi et al. [20]. Specifically, We take into account their scheme’s



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 12 of 26

special case, denoted as CM-2, which is formulated as follows:

y j = k j −
g(k j)

g′(k j)
,

z j = y j − (−6t3
j +5t2

j +2t j +1)
g(y j)

g′(k j)
,

w j = z j − (1+2t j +4u j +6t2
j + k j)

g(z j)

g′(k j)
,

k j+1 = w j − (I(t j)+ J(k j)+K(u j)+L(t j,u j) (40)

+M(p j,q j,r j)+N(t j,k j,u j,r j)).
g(w j)

g′(k j)
,

where

I(t j) = 6t2
j +2t j, J(k j) =−k

3

j + k j +1,

K(u j) = 4u j −4u2
j ,

L(t j,u j) = t ju j +6t2
j u j +2t3

j u j −10t ju2
j ,

M(p j,q j,r j) = r j +2q j +8p j,

N(t j,k j,u j,r j)) = 2t jr j +2k ju j +6t2
j r j −4k2

j u j +24t4
j u j,

t j =
g(y j)

g(k j)
, u j =

g(z j)

g(k j)
, p j =

g(w j)

g(k j)
, k j =

g(z j)

g(y j)
,

q j =
g(w j)

g(y j)
, r j =

g(w j)

g(z j)
.

Sivakumar et al. [28] introduced the following four-point iterative method utilizing the divided
difference technique, known as CM-3:

y j = k j −u(k j) ,z j = k j +g(k j)
4,

z j = k j −u(k j)[
g(k j)−g(y j)

g(k j)−2g(y j)
],

w j = z j −
g(z j)

q′
(z j)

,

k j+1 = w j −
g(w j)

r′(w j)
, (41)

where

u(k j) =
g(k j)

g′(k j)
,

q′(z j) = a1+2a2(z j − k j)+3a3(z j − k j)
2,

r′(w j) = b1 +2b2(w j − k j)+3b3(w j − k j)
2 +4b4(w j − k j)

3,

a1 = g′(k j) = b1,



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 13 of 26

a2 =
g[y j,k j,k j](z j − k j)−g[z j,k j,k j](y j − k j)

z j − y j
,

a3 =
g[z j,k j,k j]−g[y j,k j,k j]

z j − y j
,

b4 =
(g[y j,k j,k j](S3 −S2)+g[z j,k j,k j](S3 +S1)+g[w j,k j,k j](S2 −S1))

−S2
1S2 +S1S2

2 +S2
1S3 −S2

2S3 −S1S2
3 +S2S2

3
,

b3 =
(g[y j,k j,k j](S2

2 −S2
3)+g[z j,k j,k j](S2

3 −S2
1)+g[w j,k j,k j](S2

1 −S2
2))

−S2S2
1 +S2

2S1 −S3S2
2 −S2

3S1 +S3S2
1 +S2

3S2
,

b2 =

(g[y j,k j,k j](−S2
2S3 +S2S2

3)+g[z j,k j,k j](−S2
1S3 +S1S2

3)
+g[w j,k j,k j](−S2

1S2 +S1S2
2))

−S2
1S2 +S1S2

2 +S2
1S3 −S2

2S3 −S1S2
3 +S2S2

3
,

S1 = y j − k j, S2 = z j − k j, S3 = w j − k j.

We have taken another method developed by Nusrat et al. [29] for comparison, which con-
structs optimal sixteenth-order derivative-free root-finding methods based on Hermite interpola-
tion, referred to as CM-4.

y j = k j −
g(k j)

g[w j,k j]
,

w j = k j +g(k j)
4, j ≥ 0,

z j = y j −
g(y j)

2g[y j,k j]−g[w j,k j]
,

t j = z j −
g(z j)

K′
3(z j)

,

k j+1 = t j −
g(t j)

K′
4(t j)

,

where

K′
3(z) = g[z,k]

(
2+

z− k
z− y

)
− (z− k)2

(y− k)(z− y)
g[y,k]+g[w,k]

z− y
y− k

,

K′
4(t) = g[t,z]+ (t − z)g[t,z,y]+ (t − z)(t − y)g[t,z,y,k]

+(t − z)(t − y)(t − k)g[t,z,y,k,2],

g[t,z,y,k,2] =
1

(t − k)2(t − y)

(
g[t,z]−g[z,y]

)
− 1

(t − k)2(z− k)

(
g[z,y]−g[y,k]

)
− 1
(t − k)(z− k)2

(
g[y,k]−g[w,k]

)
+

1
(t − k)(z− k)(y− k)

(
g[y,k]−g[w,k]

)
.

We adopt the approach proposed by Soleymani et al. [30], which constructs an optimal 16th-
order iterative class for approximating simple zeros of nonlinear equations, referred to as CM-5 is
given as:

y j = k j −
g(k j)

g′(k j)
,



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 14 of 26

z j = y j −
g(y j)

g(k j)−2g(y j)

g(k j)

g′(k j)
, ,

w j = z j −
[

1+
g(z j)

g(k j)

]
g[k j,y j]g(z j)

g[k j,z j]g[y j,z j]
,

k j+1 = w j −
(1+b5(w j − k j))

2g(w j)

g′(k j)+2b3(w j − k j)+(3b4 +b3b5)(w j − k j)2 +2b4b5(w j − k j)3 , (42)

where

b5 =
g[z,k]WY (w− y)+Z(g[w,k](y− z)Y +g[y,k](z−w)W )− (w− y)(w− z)(y− z)g′(k)

(Y Z(z− y)g(w)+(w− z)((w− y)(y− z)g(k)+WZg(y))+WY (y−w)g(z))
,

b2 = g′(k)+g(k)b5,

b4 =
(z− k)g[y,k]+ (k− y)g[z,k]+ (y− z)b2 +((z− k)g(y)+(k− y)g(z))b5

(k− z)(k− y)(y− z)
,

b3 = g[w,z,k]+g[w,k]b5 − (w− k)b4, (43)

wherein y− k = Y , z− k = Z, and w− k =W .
Now, our focus shifts to assessing the effectiveness of the proposed schemes through their

application to various test functions and real-world problems provided below:

Example 1. Perfectly Mixed Reactor
Consider a problem involving a Perfectly Mixed Reactor [31]. In this reactor type, the concen-

tration within the tank equals the output concentration of the reactor. The objective is to determine
the concentration of a chemical within a completely mixed reactor. The objective here is to as-
certain the concentration of a chemical within a completely mixed reactor. This is governed by a
differential equation:

dc
dk

= (cin − c)
Q f

VR
. (44)

Here the inflow chemical concentration is represented by cin, while Q f denotes the flow rate of
fluid entering the reactor with a volume VR. The solution to equation (44) is given by:

c(k) = (1− e−
k
τ )cin, (45)

where τ = VR
Q f
, represents the mean residence time. Now, let’s consider the equation representing

the concentration of chemical within completely mixed reactor:

g1(k) = (1− e-0.04k)cin + coe−0.04k. (46)

Here, cin = 10 represents the inflow concentration, and co = 4 is the initial concentration of the
chemical.

The exact root of g1(k) is γ =−12.770059497670 · · · . We choose the initial approximation for
g1(k) is k∗ =−13.51. Table 1 explains the numerical results for test function g1(k). Table 1 shows
that the recently developed approaches SFHM-1 and SFHM-2 outperform than the previously
published methods when compared in terms of COC and successive iterations. These methods



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 15 of 26

Table 1: Convergence behavior for g1(k)

g1(k) = 10−6e−0.04k, k∗ =−13.51
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 7.39×10−1 1.63×10−28 3.53×10−471 16.00
SFHM-2 7.39×10−1 9.23×10−29 1.97×10−475 16.00

CM-1 7.39×10−1 7.34×10−28 7.27×10−460 15.99
CM-2 7.39×10−1 5.43×10−25 1.41×10−410 15.07
CM-3 7.39×10−1 1.00×10−29 D D
CM-4 7.08×10−15 7.08×10−15 7.08×10−15 0
CM-5 7.08×10−15 7.08×10−15 7.08×10−15 0

exhibit a computational order of convergence of 16, indicating highly accurate results. On the
other hand, the methods CM-1 and CM-2 show slightly lower performance, with computational
orders of convergence of 15.99 and 15.07, respectively. The scheme CM-3 fails to converge, as
indicated by the divergence observed in its results. Meanwhile, CM-4 and CM-5 both fail to
progress using the suggested initial values, resulting in a computational order of convergence of
zero.

Example 2. Kepler’s Equation of Celestial Motion
In the field of Celestial mechanics, Kepler’s equation [32] serves as a fundamental corner-

stone, establishing a vital link between the temporal dynamics and spatial configurations of celes-
tial bodies, notably planets, concerning a designated starting point. This equation, known for its
significance, intricately relates the mean anomaly, eccentric anomaly, and the eccentricity of the
orbital path formulated as:

g2(k) = k− esin(k)−M. (47)

where k symbolizes the eccentric anomaly, M embodies the mean anomaly and e represents the
eccentricity of the orbit. In practical astronomical scenarios, consider a satellite orbiting along
an elliptical trajectory characterized by an eccentricity of e = 0.9995 and a mean anomaly M =
0.01.Compute the eccentric anomaly while maintaining the criteria where 0 ≤ e ≤ 1 and 0 ≤ M ≤
π.

The exact root of function g2(k) is 0.3899777749463 · · · . Table 2 provides a thorough compar-
ison between the established counterparts and the suggested four-point iterative schemes SFHM-1
and SFHM-2 in order to provide more understanding of the computational approaches used. In
particular, we investigate how the initial approximation k∗ = 0.381 influence the convergence be-
havior and computational efficiency of these iterative techniques. Among the previously published
methods, CM-1 and CM-2 also exhibit a strong performance with a COC of 16.0. However, CM-3
shows a slightly reduced COC of 15.8, indicating a minor drop in accuracy. On the other hand,
CM-4 and CM-5 display significantly lower COC values of 5.53 and 3.94, respectively, highlight-
ing their comparatively weaker performance for this specific example. These results confirm the
superior efficiency and robustness of SFHM-1 and SFHM-2 in handling this nonlinear equation,
further emphasizing the importance of advanced iterative methods in this domain.

Example 3. Population Dynamics in Predator-Prey Ecosystems



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 16 of 26

Table 2: Convergence behavior for g2(k)

g2(k) = k−0.9995sin(k)−0.01, k∗ = 0.381
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 8.98×10−3 2.25×10−28 6.02×10−438 16.0
SFHM-2 0.8×10−2 2.25×10−28 5.95×10−438 16.0

CM-1 0.8×10−2 2.99×10−26 5.07×10−402 16.0
CM-2 0.8×10−2 3.86×10−23 4.30×10−350 16.0
CM-3 0.8×10−2 4.57×10−28 7.61×10−433 15.8
CM-4 5.59×10−1 5.84×10−2 6.21×10−14 5.53
CM-5 7.09×10−1 1.62×10−2 6.21×10−14 3.94

Consider a predator-prey ecosystem [33] where the prey population k, experiences both natu-
ral growth and nonlinear mortality due to overpopulation and competition. This ecosystem could
apply to various animal populations where the death rate increases nonlinearly due to overcrowd-
ing effects, resource limitations, or intensified predation as population density grows.

In such ecosystems, the growth rate of the prey population can be modeled as:

g3(k) =−20k5 − 1
2

k+
1
2
. (48)

From the five roots of the function g3(k), we choose γ = 0.427677296931003 · · · . For the
function g3(k), the initial approximation is k∗ = 0.5. Table 3 illustrates the numerical results.

Table 3: Convergence analysis for g3(k)

g3(k) =−20k5 − 1
2 k+ 1

2 , k∗ = 0.5
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 7.23×10−2 4.57×10−11 2.03×10−157 15.90
SFHM-2 7.23×10−2 1.45×10−12 3.48×10−182 15.85

CM-1 7.23×10−2 2.42×10−10 7.89×10−145 15.86
CM-2 7.23×10−2 2.57×10−8 2.78×10−109 15.65
CM-3 7.23×10−2 1.42×10−11 2.32×10−166 15.94
CM-4 1.23×10−15 1.21×10−230 0 D
CM-5 5.16×10−14 5.41×10−202 0 D

Table 3 presents the results for g3(k) with k∗ = 0.5. The proposed methods, SFHM-1 and
SFHM-2 demonstrating strong and reliable performance. Among the previously proposed meth-
ods, CM-1 and CM-3 also deliver comparable results with COC values But slower in terms of
error reduction. However, CM-2 exhibits a slightly lower COC of 15.09, indicating reduced effi-
ciency in this case. The methods CM-4 and CM-5 fail to converge, as shown by their divergence
behavior. This highlights their inability to handle the nonlinearity of this example effectively. The
results reaffirm the superior performance and reliability of the proposed methods, SFHM-1 and
SFHM-2, for this nonlinear equation while also pointing out limitations in some of the compared
schemes.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 17 of 26

Example 4. We now consider a different nonlinear test function to be:

g4(k) = (k+2)ek −1. (49)

The nonlinear function g4(k) exhibits an exact root of γ = −0.4428544010023 · · · . and k∗ =
−0.41 is our first approximation. Table 4 displays the comparison findings of g4(k). In terms

Table 4: Convergence analysis for g4(k)

g4(k) = (k+2)ek −1, k∗ =−0.41
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 3.28×10−2 2.73×10−27 1.69×10−428 15.99
SFHM-2 3.28×10−2 7.54×10−28 5.94×10−438 15.99

CM-1 3.28×10−2 2.56×10−26 5.73×10−412 15.99
CM-2 3.28×10−2 4.87×10−24 2.70×10−372 15.95
CM-3 3.28×10−2 4.05×10−28 D D
CM-4 3.02×10−26 4.26×10−411 2.00×10−2000 4.12
CM-5 9.15×10−25 1.62×10−385 1.00×10−2000 4.47

of computational order of convergence and the error between two consecutive iterations, Table
4 demonstrates that the suggested schemes SFHM-1 and SFHM-2 exhibit superior performance
compared to the other approaches presented for comparison. In contrast, the methods CM-1,
CM-2, and CM-3 show relatively weaker outcomes, with CM-3 even displaying divergence. Fur-
thermore, while CM-4 and CM-5 achieve convergence, their computational orders of convergence
remain lower than those of the proposed schemes.

Example 5. Neural Activation in Response to Complex Stimuli
The nonlinear equation:

g5(k) = sin(2cos(k))−1− k2 + esink3
. (50)

represents the activation response of a neural circuit to varying intensities or complexities of
external stimuli, where k represents the stimulus intensity or frequency. In neuroscience, neurons
or neural networks often respond to external stimuli in a complex, nonlinear fashion. Sensory
processing in the brain, such as auditory or visual inputs, may elicit responses influenced by
multiple factors, including the frequency, amplitude, and complexity of the stimulus [34].

This model can be applied to studying sensory processing, especially in response to oscillatory
or periodic stimuli, where nonlinearities play a critical role in signal processing and perception.

We choose γ =−0.7848959 · · · as the root for the function g5(k) , and the appropriate starting
approximation is k∗ = −0.7. Table 5 lists the comparative results for test function g5(k). Table 5
demonstrates that the suggested techniques SFHM-1 and SFHM-2 are not only convergent but
also more efficient computationally than the existing approaches CM-21-CM-3 in terms of error
of successive iteration for the test function g5(k). Also, the methods CM-21 and CM-3 exhibit zero
computational convergence order.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 18 of 26

Table 5: Convergence analysis for g5(k)

g5(k) = sin(2cos(k))−1− k2 + esin(k3), k∗ =−0.7
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 8.48×10−2 3.04×10−19 5.96×10−297 15.91
SFHM-2 8.48×10−2 1.80×10−16 5.11×10−252 16.05

CM-1 8.48×10−2 6.51×10−17 1.54×10−259 16.00
CM-2 8.48×10−2 2.09×10−13 6.21×10−205 16.01
CM-3 8.48×10−2 5.13×10−19 6.98×10−295 16.02
CM-4 1.63×10−15 1.62×10−51 1.62×10−51 0
CM-5 1.63×10−15 1.62×10−51 1.62×10−51 0

Example 6. Population Growth in Ecosystems with Periodic Resources
We consider the nonlinear test function to represent the growth of a biological population

where k represents a variable such as time, resource availability, or environmental conditions.

g6(k) = ln(1+ k2)+ ek sin(k). (51)

In ecology, population dynamics can be influenced by both intrinsic factors (like growth rates) and
extrinsic factors (such as resource availability, seasonal changes, or environmental conditions)
[35]. This model can be applied to studying ecosystem populations where resource availability
and environmental cycles (like seasonal changes) significantly impact growth.

The exact root of the above equation is γ = 0.0 Numerical results are shown in Table 6 with
the starting approximation k∗ = 0.01 of function g6(k). Table 6 illustrates the convergence per-

Table 6: Convergence analysis for g6(k)

g6(k) = ln(1+ k2)+ ek sin(k), k∗ = 0.01
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 9.99×10−3 1.59×10−29 4.48×10−459 16.02
SFHM-2 9.99×10−3 3.21×10−27 6.19×10−419 15.99

CM-1 9.99×10−3 7.45×10−27 9.80×10−413 15.99
CM-2 1.23×10−2 2.05×10−26 6.89×10−405 15.97
CM-3 9.99×10−3 1.51×10−28 D D
CM-4 1.26×10−28 7.83×10−443 0 D
CM-5 4.68×10−26 4.09×10−399 4.94×10−3991 9.63

formance of newly proposed methods (SFHM-1 and SFHM-2) compared to established methods
(CM-1, CM-2, CM-3, CM-4, and CM-5) for the nonlinear equation g6(k). The proposed methods,
SFHM-1 and SFHM-2, demonstrate remarkable performance with significantly smaller error
levels and higher COC values of 16.02 and 15.99, respectively, outperforming all other methods.
Among the established methods, CM-1 and CM-2 also achieve required COC values but their er-
ror levels are slightly higher than those of the proposed techniques. The methods CM-3 and CM-4
encounter divergence and CM-5 exhibits the weakest performance with a low COC of 9.63 and
larger error levels.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 19 of 26

Example 7. We use the following nonlinear test function:

g7(k) = (k−1)3 −1. (52)

With three approximate roots in the preceded equation, we choose γ = 2 We chose k∗ = 2.1 as
our first approximation. The numerical findings are displayed in Table 7. When comparing con-

Table 7: Convergence analysis for g7(k)

g7(k) = (k−1)3 −1, k∗ = 2.1
Methods |k2 − k1| |k3 − k2| |k4 − k3| COC
SFHM-1 9.99×10−2 7.46×10−18 1.27×10−275 15.98
SFHM-2 9.99×10−2 3.75×10−17 1.92×10−263 15.96

CM-1 9.99×10−2 3.13×10−16 1.03×10−247 15.96
CM-2 9.99×10−2 1.13×10−16 1.34×10−252 15.78
CM-3 9.99×10−2 7.85×10−18 4.16×10−275 15.97
CM-4 6.38×10−14 1.83×10−210 0 D
CM-5 2.94×10−15 5.66×10−231 0 D

cerning computational order of convergence and absolute error, the suggested schemes SFHM-1
and SFHM-2 outperform CM-1-CM-5, as demonstrated by Table 7 for the function g7(k).

5. Dynamical Investigation of Proposed Family

Basins of attraction are used to study the dynamical behavior of a rational function connected
to an iterative approach, which offers essential perspectives on the convergence and stability of the
method. It was initially proposed by Vrscay and Gilbert [36] to analyze the complicated dynamics
of iteration systems. Later, several researchers used this strategy in their writings [37, 38]. A
graphic that illustrates the behavior of such an algorithm as a function of the different beginning
positions is called the basin of attraction.

Let a rational mapping represented by χ : C→ C on the complex plane, where C represents
the Riemann sphere. The orbit of the point ρo ∈ C is represented by the set that follows as,{

ρo,χ(ρo),χ
2(ρo), . . . ,χ

n(ρo), . . .
}
.

If a point ρ0 ∈C fulfills χn(ρ0) = ρ0, where n is the smallest value that satisfies this condition,
then the point is termed a periodic point with a minimal period n. A periodic point with a minimum
period of one is called a fixed point of χ . Fixed points are categorized based on the multiplier
|χ ′(ρ0)|. Taking into consideration the associated multiplier, a fixed point ρ0 can be classified as
follows:

• If |χ ′(ρ0)|< 1, then the point is an attractor.

• If |χ ′(ρ0)|> 1, then the point is repelling.

• If |χ ′(ρ0)|= 1, then the point is neutral.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 20 of 26

• If |χ ′(ρ0)|= 0, then the point is a superattractor.

The collection of preimages of any order is referred to as C(α) if α is a rational function that
attracts fixed points.

R(α) = {ρ0 ∈ C | χ
n(ρ0)→ α,n → ∞} .

A set of points known as the Fatou set is where orbits converge while attempting to attract a
fixed point. The Julia set, the closure of a collection of repelling fixed points is its counterpart and
establishes the limits separating the attraction basins.

(a) SFHM-1 (b) SFHM-1 (c) SFHM-2 (d) SFHM-2

(e) CM-1 (f) CM-1 (g) CM-2 (h) CM-2

(i) CM-3 (j) CM-3 (k) CM-4 (l) CM-4

(m) CM-5 (n) CM-5

Figure 1: Dynamical behavior of different methods for z3 (ρ)

To establish a basin of attraction, we take two distinct approaches. The first approach utilizes
a rectangular box defined by [−2,2]× [−2,2] ∈ C. Iterations are performed up to a maximum



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 21 of 26

(a) SFHM-1 (b) SFHM-1 (c) SFHM-2 (d) SFHM-2

(e) CM-1 (f) CM-1 (g) CM-2 (h) CM-2

(i) CM-3 (j) CM-3 (k) CM-4 (l) CM-4

(m) CM-5 (n) CM-5

Figure 2: Dynamical behavior of different methods for z4 (ρ)

of 25, and the root of the sequence generated using the iterative technique converges with the
stopping condition |g(u j)| ≤ 10−5. In this approach, the points are darkened based on the number
of iterations to visualize the convergence speed. Depending on the root of convergence of the
iterative algorithm, we assign a color to each initial guess ρ0. Dark colors are assigned to indicate
divergence.

In the second approach, each root is given a distinct shade, and the points are also colored
based on the distance to the closest root. In other words, after identifying the nearest root, we
assign its color to the starting point. Consequently, this approach illustrates which of the basic
iterative methods converges. We use an error estimate of less than 10−5 and a maximum of 25 it-
erations. With this method, each starting estimate is assigned a unique shade, and the convergence
of finding a root of a specific nonlinear polynomial determines the number of iterations.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 22 of 26

(a) SFHM-1 (b) SFHM-1 (c) SFHM-2 (d) SFHM-2

(e) CM-1 (f) CM-1 (g) CM-2 (h) CM-2

(i) CM-3 (j) CM-3 (k) CM-4 (l) CM-4

(m) CM-5 (n) CM-5

Figure 3: Dynamical behavior of different methods for z5 (ρ)

Two distinct strategies were utilized to create polynomiographs using MATLAB R2014a. To
create dynamic planes, we consider three complex functions for dynamic testing. These functions
are provided by:

z3 (ρ) = ρ
3 −1, ρ = 1.0, −0.5000+0.86605I, −0.5000−0.86605I,

z4 (ρ) = ρ
4 −10ρ

2 +9, ρ = 1.0, −3.0, 3.0, −1.0,

z5(ρ) = ρ
5 −1, ρ = 1.0, 0.3090−0.951I, 0.309+0.951I,

−0.809−0.587I, 0.809+0.587I

The complex planes of the suggested techniques (SFHM-1) and (SFHM-2) as well as the previ-
ously constructed algorithms by Sharifi et al. [20] (CM-1), Sharma et al. [21] (CM-2), Sivakumar



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 23 of 26

et al. [28] (CM-3), Nusrat et al. [29] (CM-4), and Soleymani et al. [30] (CM-5) are shown in
Figures 1-3. There are two different kinds of basins of attraction illustrated in all these figures:

Based on the larger and darker regions of convergence observed in Figures 1-3 when compared
to other techniques, it is evident that our proposed methods, SFHM-1 and SFHM-2, exhibit a
significantly higher order of convergence. These methods also encompass a broader range of
initial conditions that successfully lead to convergence. Furthermore, the convergence regions of
SFHM-1 and SFHM-2 are not only wider but also display substantially reduced chaotic behavior
compared to the competing methods CM-1, CM-2, and CM-3. The methods CM-4 and CM-5
exhibit extensive divergence regions across all tested complex functions, indicating their limited
reliability and effectiveness in these scenarios.

For the complex polynomial z3(ρ), SFHM-1 and SFHM-2 exhibit larger and darker conver-
gence regions with significantly reduced chaotic behavior compared to CM-2, CM-4, and par-
ticularly CM-5, which demonstrates an extensive divergence region. Although CM-1 and CM-3
display relatively smaller divergence regions, they exhibit higher chaotic behavior, making them
less stable in this case. A similar pattern is observed for the complex polynomial z4(ρ), where
SFHM-1 and SFHM-2 consistently demonstrate superior stability and broader convergence re-
gions. In the case of z5(ρ), these methods continue to perform effectively, maintaining a large
convergence region with minimal chaotic behavior and reduced divergence. In contrast, the other
methods exhibit significantly larger divergence regions and pronounced chaotic behavior, leading
to slower convergence and instability in initial conditions.

Consequently, the results indicate that SFHM-1 and SFHM-2 maintain their stability even as
the degree of the complex polynomial increases, confirming their effectiveness for high-degree
polynomials. Their ability to provide larger, more stable convergence regions while ensuring
accuracy and efficiency establishes them as superior alternatives to the existing techniques CM-1-
CM-5.

6. Conclusion

In this research, we have presented novel four-point iterative methods, referred to as SFHM-
1 and SFHM-2, which are developed within the framework of a weighted Steffensen-Ostrowski
family. These methods utilize the Lagrange interpolation technique, involving five function eval-
uations per iteration, achieving a validated sixteenth-order convergence and supporting the Kung-
Traub conjecture. Additionally, we have demonstrated the practical advantages of these methods
by applying them to complex, real-world scenarios, including Kepler’s celestial motion, ideally
mixed reactors, and models of predator-prey dynamics, neural activation, and periodic ecosystem
growth. Compared to existing iterative techniques such as CM-1, CM-2, CM-3, CM-4,and CM-5,
our methods exhibit superior convergence speed and efficiency. A detailed basin of attraction anal-
ysis further illustrates the robustness and adaptability of SFHM-1 and SFHM-2 across a diverse
range of scenarios. The consistency between our theoretical predictions and numerical results un-
derscores the practical relevance and reliability of the proposed methods. Our findings suggest that
the SFHM-1 and SFHM-2 schemes provide a competitive alternative for solving nonlinear equa-
tions, with potential applications across various complex problem domains, encouraging further
research and practical exploration.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 24 of 26

Future directions in this field involve extending higher-order iterative methods to complex,
multi-dimensional problems and systems of equations. Exploring hybrid approaches that integrate
machine learning for parameter estimation holds significant potential for boosting convergence
rates and enhancing robustness. Additionally, further theoretical investigations into stability, error
dynamics, and adaptation to diverse nonlinear problems will contribute to the broader applicabil-
ity of these methods. This progress could lead to more efficient and reliable numerical solvers
capable of addressing the growing complexity of real-world applications across various scientific
and engineering disciplines.

Acknowledgements

Fouzia Amir expresses sincere gratitude to Asia International University, Bukhara, Uzbekistan
for its financial support. Muhammad Bilal Riaz acknowledges the financial support of the Euro-
pean Union under the REFRESH – Research Excellence For Region Sustainability and High-tech
Industries project number CZ.10.03.01/00/22 003/0000048 via the Operational Programme Just
Transition.

References

[1] A. M. Ostowski. Solution of Equations and System of Equations. Academic Press, 1960.
[2] H. T. Kung and J. F. Traub. Optimal order of one-point and multipoint iteration. Journal of

the Association for Computing Machinery, 21(4):643–651, 1974.
[3] M. S. Petković, B. Neta, L. D. Petković, and J. D. Džunić. Multipoint Methods for Solving

Nonlinear Equations. Elsevier, 2013.
[4] J. F. Traub. Iterative Methods for the Solution of Equations. Prentice-Hall, 1964.
[5] W. Bi, H. Ren, and Q. Wu. Three-step iterative methods with eight-order convergence

for solving nonlinear equations. Journal of Computational and Applied Mathematics,
225(1):105–112, 2009.

[6] D. K. R. Babajee, K. Madhu, and J. Jayaraman. A family of higher-order multi-point iterative
methods based on power mean for solving nonlinear equations. Afrika Matematika, 27(5–
6):865–876, 2016.

[7] Y. Chu, N. Rafiq, M. Shams, S. Akram, N. A. Mir, and H. Kalsoom. Computer methodolo-
gies for the comparison of some efficient derivative-free simultaneous iterative methods for
finding roots of nonlinear equations. Computers, Materials and Continua, 66(1):275–290,
2020.

[8] G. Thangkhenpau, S. Panday, L. C. Bolunduţ, and L. Jäntschi. Efficient families of multi-
point iterative methods and their self-acceleration with memory for solving nonlinear equa-
tions. Symmetry, 15(8):Article 1546, 2023.

[9] S. Panday, A. Sharma, and G. Thangkhenpau. Optimal fourth and eighth-order itera-
tive methods for nonlinear equations. Journal of Applied Mathematics and Computing,
69(1):953–971, 2023.

[10] M. Q. Khirallah and A. M. Alkhomsan. A new fifth-order iterative method for solving non-



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 25 of 26

linear equations using the weight function technique and the basins of attraction. Journal of
Mathematics and Computer Science, 28(3):281–293, 2023.

[11] S. Abdullah, N. Choubey, and S. Dara. An efficient two-point iterative method with mem-
ory for solving nonlinear equations and its dynamics. Journal of Applied Mathematics and
Computing, 70(1):285–315, 2024.

[12] M. U. D. Junjua, F. Zafar, and N. Yasmin. Optimal derivative-free root finding methods
based on inverse interpolation. Mathematics, 7(2):164, 2019.

[13] M. U. D. Junjua, S. Abdullah, M. Kansal, and S. Ahmad. On traub–steffensen-type iteration
schemes with and without memory: Fractal analysis using basins of attraction. Fractal and
Fractional, 8(12):698, 2024.

[14] S. Abdullah, N. Choubey, S. Dara, M. U. D. Junjua, and T. Abdullah. A robust and op-
timal iterative algorithm employing a weight function for solving nonlinear equations with
dynamics and applications. Axioms, 13(10):675, 2024.

[15] S. Abdullah, N. Choubey, and S. Dara. Optimal fourth-and eighth-order iterative methods
for solving nonlinear equations with basins of attraction. Journal of Applied Mathematics
and Computing, pages 1–31, 2024.

[16] S. Akram, M. Khalid, M. u. D. Junjua, S. Altaf, and S. Kumar. Extension of king’s iterative
scheme by means of memory for nonlinear equations. Symmetry, 15(5):1116, 2023.

[17] S. Abdullah, N. Choubey, and S. Dara. Two novel with and without memory multi-point
iterative methods for solving nonlinear equations. Communications in Mathematics and
Applications, 15(1):9, 2024.

[18] D. K. R. Babajee and R. Thukral. On a four-point sixteen-order king family of iterative
method for solving nonlinear equations. International Journal of Mathematics and Mathe-
matical Sciences, 2012:Article ID 979245, 13 pages, 2012.

[19] F. Zafar, N. Hussain, Z. Fatimah, and A. Kharal. Optimal sixteenth-order convergent method
based on quasi-hermite interpolation for computing roots. The Scientific World Journal,
2014:Article ID 930941, 7 pages, 2014.

[20] S. Sharifi, M. Salimi, S. Siegmund, and T. Lotfi. A new class of optimal four-point methods
with convergence order 16 for solving nonlinear equations. Mathematics and Computers in
Simulation, 119:69–90, 2016.

[21] J. R. Sharma and S. Kumar. Efficient methods of optimal eighth and sixteenth-order conver-
gence for solving nonlinear equations. SeMA Journal, 75(2):229–253, 2018.

[22] R. Behl, S. Amat, Á. A. Magreñán, and S. S. Motsa. An efficient optimal family of sixteenth-
order methods for nonlinear models. Journal of Computational and Applied Mathematics,
354:271–285, 2019.

[23] J. R. Sharma and H. Arora. Efficient ostrowski-like methods of optimal eighth- and sixteenth-
order convergence and their dynamics. Afrika Matematika, 30(5–6):921–941, 2019.

[24] D. Ćebić, N. M. Ralević, and M. Marčeta. An optimal sixteenth-order family of methods for
solving nonlinear equations and their basins of attraction. Mathematical Communications,
25(2):269–288, 2020.

[25] S. Akram, H. Khalid, T. Rasulov, M. Khalid, and M. U. Rehman. A novel variation of
16th-order iterative scheme for models in blood stream, chemical reactor, and its dynamics.
Lobachevskii Journal of Mathematics, 44(9):5116–5131, 2023.



S. Akram et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5822 26 of 26

[26] V. Kanwar, R. Bala, and M. Kansal. Some new weighted eighth-order variants of steffensen-
king’s type family for solving nonlinear equations and its dynamics. SeMA Journal,
74(1):75–90, 2017.

[27] M. U. D. Junjua, S. Akram, and A. Tariq. Construction of optimal derivative-free iterative
methods for nonlinear equations using lagrange interpolation. Journal of Prime Research in
Mathematics, 16(1):30–45, 2020.

[28] P. Sivakumar, K. Madhu, and J. Jayaraman. Optimal eighth and sixteenth-order iterative
methods for solving nonlinear equation with basins of attraction. Applied Mathematics E-
Notes, 21:320–343, 2021.

[29] N. Yasmin, F. Zafar, and S. Akram. Optimal derivative-free root finding methods based on
the hermite interpolation. Journal of Nonlinear Sciences and Applications, 9:4427–4435,
2016.

[30] F. Soleymani, S. Shateyi, and H. Salmani. Computing simple roots by an optimal sixteenth-
order class. Journal of Applied Mathematics, 2012(1):958020, 2012.

[31] M. S. Petković, B. Neta, L. D. Petković, and J. Džunić. Multipoint methods for solving
nonlinear equations: A survey. Applied Mathematics and Computation, 226:635–660, 2014.

[32] A. Cordero, J. L. Hueso, E. Martı́nez, and J. R. Torregrosa. Generating optimal derivative-
free iterative methods for nonlinear equations by using polynomial interpolation. Mathemat-
ical and Computer Modelling, 57(7-8):1950–1956, 2013.

[33] S. H. Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chem-
istry, and engineering. CRC Press, 2018.

[34] E. R. Kandel, J. H. Schwartz, T. M. Jessell, S. Siegelbaum, A. J. Hudspeth, and S. Mack.
Principles of neural science, volume 4. McGraw-Hill, New York, 2000.

[35] W. E. Grant and T. M. Swannack. Ecological modeling: a common-sense approach to theory
and practice. John Wiley & Sons, 2011.

[36] E. R. Vrscay and W. J. Gilbert. Extraneous fixed points, basin boundaries and chaotic dynam-
ics for schröder and könig rational iteration functions. Numerische Mathematik, 52(1):1–16,
1988.

[37] J. L. Varona. Graphic and numerical comparison between iterative methods. The Mathemat-
ical Intelligencer, 24(1):37–46, 2002.

[38] A. Naseem, M. A. Rehman, and T. Abdeljawad. Computational methods for nonlinear equa-
tions with some real-world applications and their graphical analysis. Intelligent Automation
and Soft Computing, 30(3):545–558, 2021.