A New Three-Step Root-Finding Numerical Method and Its Fractal Global Behavior

: There is an increasing demand for numerical methods to obtain accurate approximate solutions for nonlinear models based upon polynomials and transcendental equations under both single and multivariate variables. Keeping in mind the high demand within the scientiﬁc literature, we attempt to devise a new nonlinear three-step method with tenth-order convergence while using six functional evaluations (three functions and three ﬁrst-order derivatives) per iteration. The method has an efﬁciency index of about 1.4678, which is higher than most optimal methods. Convergence analysis for single and systems of nonlinear equations is also carried out. The same is veriﬁed with the approximated computational order of convergence in the absence of an exact solution. To observe the global fractal behavior of the proposed method, different types of complex functions are considered under basins of attraction. When compared with various well-known methods, it is observed that the proposed method achieves prespeciﬁed tolerance in the minimum number of iterations while assuming different initial guesses. Nonlinear models include those employed in science and engineering, including chemical, electrical, biochemical, geometrical, and meteorological models.


Introduction
The study of iterative methods for solving nonlinear equations and systems appears to be a very important area in theory and practice. Such problems appear not only in applied mathematics but also in many branches of science including engineering (design of an electric circuit), physics (pipe friction), chemistry (greenhouse gases and rainwater), biology (steady-state of Lotka-Volterra system), fluid dynamics (combined conductiveradiative heat transfer), environmental engineering (oxygen level in a river downstream from a sewage discharge), finance (option pricing), and many more. The study of nonlinear models is a vital area of research in numerical analysis. Interesting applications in pure and applied sciences can be studied in the general construction of the nonlinear equations expressed in the form g(x) = 0. Due to their significance, several iterative methods have been devised under certain situations since it is near to impossible to obtain exact solutions of models that are nonlinear in nature. These iterative methods have been constructed using different existing methods such as Taylor expansion, the perturbation method, the homotopy perturbation method, Adomian decomposition method, quadrature formula, multi-point iterative methods, the Steffensen-type methods adapted to multidimensional cases, and the variational iteration method. For detailed information, see [1][2][3][4][5][6][7][8]. Among existing iterative methods, the optimal methods are considered those that satisfy the condition for an order of convergence of 2 k−1 , where k stands for the number of function evaluations per iteration as suggested in [9]. In this way, Newton's classical method x n+1 = x n − g(x n ) g (x n ) is the optimal method with quadratic convergence. Various attempts have been made to improve the efficiency of Newton's classical method in past and recent research, as can be seen in [10][11][12][13][14][15] and most of the references cited therein.

Materials and Methods
In a general form, the uni-variate nonlinear equation can be expressed as g(x) = 0, where x is the desired quantity. It is extremely difficult to solve the nonlinear equation to find the value of x. Therefore, we attempt to devise a new, highly convergent iterative method to obtain an accurate approximate x that could yield the smallest possible error in the numerical solution. Before we continue with a discussion and derivation of the proposed method, we present some of the methods that are frequently used in the available literature. Later, we use these methods to compare the results obtained under these methods and the results obtained via the method we plan to propose.

Some Existing Methods
The iterative method ,called the Newton Rahpson method NR 2 [1,16,17] with quadratic convergence, is shown below and uses two function evaluations: one for the function g(x) itself and 1 for the first derivative g (x): x n+1 = x n − g(x n ) g (x n ) , n = 0, 1, 2, . . . , (1) where g (x n ) = 0. In [2], the authors proposed an iterative method with fifth-order convergence as abbreviated by MHM 5 . The method requires four function evaluations per iteration: two for the function itself and two first derivatives. The computational steps for the two-step method MHM 5 is described below: x n+1 = y n − g (y n ) + 3g (x n ) 5g (y n ) − g (x n ) . g(y n ) g (x n ) , n = 0, 1, 2, . . . , (2) where g (x n ) = 0 and 5g (y n ) = g (x n ). An efficient three-step iterative method with sixth-order convergence is proposed in [3]. This method is the combination of two different methods from [1,18] with second and third-order convergence, respectively. The method requires five function evaluations: three evaluations of the function itself and two evaluations of the first-order derivative per iteration. We represent the method as HM 6 . The computational steps of the method can be described as follows: z n = y n − g(y n ) g (y n ) , x x+1 = y n − g(y n ) + g(z n ) g (y n ) . n = 0, 1, 2, . . . , The three-step method with eighth-order convergence as denoted by WO 8 is proposed in [19]. The method requires four function evaluations: three evaluations of the function itself and one evaluation of the first-order derivative per iteration. The computational steps of WO 8 can be described as follows: 8g(y n ) 4g(x n ) − 11g(y n ) + 1 + g(z n ) g(y n ) . n = 0, 1, 2, . . . , (4) The iterative method with ninth-order convergence can be seen in [20]. The method requires five function evaluations: three evaluations of function itself and two evaluations of the first-order derivative per iteration. This method is abbreviated as N M 9 . The computational steps of the method are described as follows: The predictor-corrector modified Householder's method with tenth order convergence, as denoted by MH 10 , is proposed in [21]. The method is free from the second derivative and requires only five function evaluations per iteration: three evaluations of the function itself and two evaluations of the first-order derivative. The computational steps of MH 10 can be described as x n+1 = z n − g(z n ) g[z n , y n ] + (z n − y n )g[z n , y n , y n ] , where

Proposed Iterative Method
There are many recent research studies wherein researchers have presented modified iterative methods to solve nonlinear models of the form g(x) = 0. In some of these methods, modification is based on the idea of combining two existing methods to develop a new method with a better order of convergence. After being motivated by such an idea as observed in [3,[22][23][24], we have developed a new method with tenth-order convergence via the blending of two different iterative methods with second (Newton method) and fifth-order (modified Halley method) convergence as given in [1,2], respectively. We recommended the higher-order convergent method of the convergence order 2 × 5 = 10. The choice of the methods in the present work is suitable because the resultant iterative method with tenth order convergence uses only 6 function evaluations (3 function evaluations and three evaluations of the first-order derivative) per iteration. It may be noted that the choice of blending of methods is extremely important to avoid extra function evaluations that could bring additional computational cost, as can be seen in [25,26]. Hence, the proposed iterative method not only confirms higher-order convergence but also employs fewer function evaluations, as can be described by the following proposed three-step method abbreviated as PM 10 : The proposed iterative three-step method given in (8) is discussed in the flowchart presented in Figure 1. E=x1-x0 Is |E| > ? Furthermore, the efficiency index e is also computed for the proposed iterative 55 method as 10 1/6 ≈ 1.4678 whereas it would generally be shown as 10 1/3(n+n 2 ) for 56 n ≥ 1. The following Table 1 and Figure 2 can be consulted with for computation and 57 comparison of all iterative methods taken for comparison in the present research work.

58
Although function evaluations (FV) of PM 10 in the Table 1 seem to be more than some of 59 the methods under consideration but to achieve the desired accuracy, the performance 60 Figure 1. Flowchart of the proposed three-step method given in (8).
Furthermore, the efficiency index e is also computed for the proposed iterative method as 10 1/6 ≈ 1.4678, whereas it would generally be shown as 10 1/3(n+n 2 ) for n ≥ 1. The following Table 1 and the Figure 2 can be consulted for the computation and comparison of all iterative methods taken for comparison in the present research work. Although the function evaluations (FV) of PM 10 in the Table 1 seem to be more than some of the Furthermore, the efficiency index e is also computed for the proposed iterative method as 10 1/6 ≈ 1.4678, whereas it would generally be shown as 10 1/3(n+n 2 ) for n ≥ 1. The following Table 1 and the Figure 2 can be consulted for the computation and comparison of all iterative methods taken for comparison in the present research work. Although the function evaluations (FV) of PM 10 in the Table 1 seem to be more than some of the methods under consideration, to achieve the desired accuracy regarding the performance of the latter under different initial guesses (IG), the number of iterations (N) and CPU time (seconds) are better than most of the methods. This discussion is presented in Section 5.

Convergence Analysis
This section has been devoted to the proof of local convergence analysis for the proposed tenth-order method under both scalar and vector (systems) cases. Single and multivariate Taylor's series expansion have been used to obtain the required order of local convergence. It is worth noting that the convergence analysis is addressed in a similar manner to many other existing articles, and the main interest in developing higher-order methods is of the academic type. Even if higher-order methods are more complicated, the efficiency can be measured, and this is why we have included the CPU time, as found in Section 5. The theorems stated below are later used for the theoretical analysis of the convergence. Theorem 1. (Single real variable Taylor's series expansion): Suppose that r ≥ 1 is an integer and further suppose that g : R → R is an r-times differentiable function at some finite point α ∈ R. Then, there exists the following expression: where R r (x) is the remainder term, whose integral form is Theorem 2. (Multivariable Taylor's series expansion): Suppose that G : P ⊆ R n → R n is an r-times Frechet differentiable system of functions in a convex set P ⊆ R n ; then, for any x and k ∈ R n the equation given below is true: where ||R r || ≤ 1 r! sup 0<t<1 ||G (r) (x + kt)||||k|| r and

Convergence under Scalar Case
In this subsection, we theoretically prove the local order of convergence for PM 10 .
Theorem 3. Suppose that α ∈ P is the required simple root for a differentiable function g : P ⊆ R → R within an open real interval P. Then, the proposed three-step numerical method (8) possesses tenth-order convergence, and the asymptotic error term is given by Proof. Suppose α is the root of g(x n ), where x n is the nth approximation for the root by the proposed method (8), and x n = x n − α is the error term in variable x at the nth iteration step. Employing the single real variable Taylor's series given in the theorem (1) for g(x n ) around α, we obtain Similarly, using the Taylor's series for 1/g (x n ) around α, we obtain Multiplying (13) and (14), we obtain Now, substituting (15) in the first step of (8), we obtain where y n = y n − α. Using the Taylor's series (1) for g(y n ) around α, we obtain g(y n ) = g (α) y n + 1 2! g (α) 2 Similarly, using the Taylor's series for 1 g (y n ) around α, we obtain Multiplying (17) and (18), we obtain Now, substituting (19) in the second step of (8), we obtain where z n = z n − α. Using the Taylor's series (1) for g(z n ) around α, we obtain Using the Taylor's series (1) for g (y n ) around α, we obtain Using the Taylor's series (1) for g (z n ) around α, we obtain Expanding the Taylor series 1 5g (z n ) − g (y n ) and using Equations (22) and (23), we obtain Finally, substituting (24) in the third step of (8) and using Equations (16), (18), (20) and (21), we obtain Hence, the tenth-order convergence of the proposed method PM 10 given by (8) for g(x) = 0 is proved.

Convergence under Vector Case
This subsection extends the proof for the tenth-order convergence of the proposed method PM 10 given in (8) regarding solving the system of nonlinear functions G(x) = 0, where G = [g 1 (x), g 2 (x), . . . , g n (x)] from R n to R n . For the system of nonlinear functions, PM 10 can be described as follows: We present the following theorem to obtain the asymptotic error term and the order of convergence for PM 10 .
Theorem 4. Let the function G : P n ⊂ R n → R n be sufficiently differentiable in a convex set P n containing the zero α of G(x). Let us consider that G (x) is continuous and nonsingular in α. Then, the solution x obtained by using proposed three-step method PM 10 converging to α has tenth-order convergence, if an initial guess x 0 is chosen close to α.
Proof. Suppose that xn = ||x n − α||. Now, using the Taylor series described in the Theorem (2) for G(x n ), we obtain Similarly, using the Taylor's series for G (x n ) around α, we obtain Employing the Taylor's series for the inverted Jacobian matrix G (x n ) −1 around α, we obtain Multiplying (27) and (29), we obtain Now, substituting (30) in the first step of (26), we obtain Similarly, employing the Taylor series for G(y n ) and G (y n ) about α, and also for the inverted Jacobian matrix G (y n ) −1 in the second step of (26), we obtain Using the Taylor series for G (z n ) and G (y n ) around α, we obtain the following: Using the above two equations and the inverted Jacobian matrix for 5G (z n ) − Finally, substituting (35) and other required expansions in the third step of (26), we obtain Hence, the tenth-order convergence of the proposed method PM 10 given by (26) for a non-linear system of equations G(x) = 0 is proved.

Computational Estimation of the Convergence Order
When a new iterative method is proposed to compute an approximate solution for g(x) = 0, in order to verify its theoretical local order of convergence, one needs to use a parameter called the Computational Order of Convergence (COC). However, this is possible only when we have information about the exact root α for g(x) = 0. Thus, the following parameters can alternatively be employed under some constraints as described below.
Approximated Computational Order of Convergence [27]: Computational Order of Convergence [28]: Extrapolated Computational Order of Convergence [29]: where Petkovic Computational Order of Convergence [30]: All of the above formulas can be used to test the convergence order, but in the present study, ACOC as given by (37) is used since the number of iterations taken by the method PM 10 (8) is at least four in various numerically tested problems, as discussed in Section 5. Additionally, this is the best-known approach employed in most of the recently conducted research studies to verify the theoretical order of convergence.

Basins of Attraction
The stability of solutions (roots) for the nonlinear function g(z) = 0 using an iterative method can be analyzed with the help of a concept called the basins of attraction. Basins of attractions are phase-planes that demonstrate iterations employed by the iterative method, which can assume different choices for the initial guess z 0 . Such 2D regions are esthetically so beautiful that their applications are not only found in applied mathematics, but also people working in fields such as architecture, arts, and design also use the concept to obtain pleasing designs. Many other fields of applications for these basins can be seen in [31][32][33][34][35][36] and most of the references cited therein. It may be noted that linear models do not depict such dynamically eye-catching behavior, whereas the non-linearity results in features such as those seen herein under the proposed iterative method PM 10 . MATLAB's built-in routines, including contour, colormap, and color bar, are utilized to obtain the basins of attraction in the present study. In this connection, a squared region R on [−4, 4] × [− 4,4] containing 2000 by 2000 mesh points is selected for the selected functions in complex form. Some of these complex-valued functions are taken as follows for the illustration of the regions by the proposed method.
To maintain diversity, different kinds of functions including polynomials and transcendental functions are used. The regions are achieved with a tolerance of = 10 −2 and a maximum number of iterations allowed of n = 12. As can be seen in Figures 3-8, for the Example (1), the maximum number of iterations is needed by initial guesses that reside near boundaries of the regions, whereas if z 0 lies within the neighborhood of the required solution, the proposed method PM 10 does not require as many iterations as those needed by most of the methods found in the literature. The average time required to produce the dynamical planes shown in Figures 3-8 is stored in the Table 2. As expected, the complex functions that have exponential and hyperbolic components are computationally expensive.

Numerical Experiments with Discussion
Various types of test problems are considered from different sources including [3,19,22]. The approximate solutions x up to 50 decimal places are shown against each test function. The error tolerance | N | = |x N − x | to stop the number of iterations is set to 10 −200 , whereas the precision is chosen to be as large as 4000, and N shows the total number of iterations taken by the method to achieve the required tolerance. In addition, physically applicable nonlinear models such as the Van der Waals equation, the Shockley ideally diode electric circuit model, the conversion of substances in a batch reactor, and the Lorenz equations in meteorology are taken into consideration to demonstrate the applicability of the proposed method PM 10 . Obtained numerical results are tabulated for further analysis. Each table contains different initial guesses, numbers of iterations required by a method to achieve the preset error tolerance, function evaluations needed for each method, For the test problem 2 (g 1 (x)), two initial guesses are chosen for simulations, as can be seen in Table 3. For the initial guess x 0 = 1.5, it is observed that the minimum number of iterations is taken by the methods PM 10 and WO 8 to achieve the error tolerance; however, the smallest error is given by PM 10 while consuming an equivalent amount of CPU time. The method NR 2 takes the maximum number of iterations at x 0 = 1.5. Under the second initial guess x 0 = 2.0, although many methods including PM 10 take the equal number of iterations to achieve = 10 −200 , the smallest absolute error and thus smallest absolute functional value is achieved by PM 10 . This shows that if an initial guess lying near to the solution of g(x) = 0 is passed to PM 10 , then the method yields the smallest error.
For the test problem 2 (g 2 (x)), two initial guesses are chosen for simulations as can be seen in Table 4. One of them is taken far away from the approximate solution of the quintic equation g 2 (x). For x 0 = −3.8, the smallest possible absolute error is yielded by WO 8 , but at the cost of the maximum number of iterations and largest amount of CPU time. Next comes HM 6 with an absolute error of 9.2097 × 10 −670 and greater time efficiency, but it requires one more iteration when compared with PM 10 , which achieves an absolute error of 1.9515 × 10 −572 with only four iterations while consuming a reasonable amount of CPU time. When an initial guess lying near to the root is chosen, the method N M 9 achieves the smallest error, but it requires one extra iteration when compared with PM 10 and WO 8 . The most expensive methods (in terms of N, FV, and CPU time) for this particular test problem seem to be WO 8 and NR 2 .
For the transcendental problem 2 (g 3 (x)), two initial guesses are chosen for simulations, as can be seen in Table 5. Under both of the initial guesses, the maximum numbers of function evaluations are taken by HM 6 followed by PM 10 to achieve the desired error tolerance. The smallest absolute functional values are obtained with HM 6 and PM 10 , where NR 2 seems to be the most expensive method in terms of the number of iterations. Although the method MH 10 consumes the fewest number of CPU seconds, its ACOC is only eight, contradicting the theoretical order of convergence found in [21].
For the transcendental problem 2 (g 4 (x)), two initial guesses are chosen for simulations as can be seen in Table 6. One of the guesses is taken far away from the approximate solution of g 4 (x). Under both initial guesses, it is observed that the method MHM 5 takes the fewest iterations, fewest function evaluations, and fewest CPU seconds with unsatisfactory absolute errors under x 0 = −9.5. The method WO 8 diverges for the second initial guess x 0 = −9.5, whereas NR 2 is the most expensive method concerning N and FV, in particular. The proposed method PM 10 performs reasonably well under the initial guesses and does not diverge under any situation.
For the test problem 2 (g 5 (x)), three initial guesses are chosen for simulations as can be seen in Table 7. Two of the guesses lie far away from the approximate root of g 5 (x). It is easy to observe that the method PM 10 performs better than other methods, even when the initial guesses are not near to the approximate root, since the number of iterations to attain the required accuracy is the smallest with PM 10 . Once again, the most expensive method regarding N and FV proves to be NR 2 , whereas WO 8 does not succeed towards the desired root when the initial guesses are assumed to be away from the root. The absolute error achieved by PM 10 with x 0 = 2.9 is the smallest when compared with the results of other methods.
After some simplifications, one obtains the following polynomial of nonlinear form: The Van der Waals equation of state was formulated in 1873, with two constants a and b (Van der Waals constants) determined from the behavior of a substance at the critical point. The equation is based on two effects, which are the molecular size and attractive force between the molecules. The model (43) is a modified version of the ideal gas equation V = RT/nP, where n shows the number of moles, R stands for the universal gas constant (0.0820578), T is the absolute temperature, V shows the volume, and P denotes the absolute pressure. If V = 1.4 moles of benzene vapor form under P = 40 atm with a = 18 and b = 0.1154, then one has The approximate solution up to 50 dp is given as: Being cubic, the Equation (45) certainly possesses one real root. Here, we aim to show the performance of PM 10 on this model. Therefore, the model is numerically solved by PM 10 and the other six methods chosen for comparison. It can be observed in Table 8 that PM 10 achieves the smallest possible error along with functional values nearest to 0 in a reasonable amount of time, irrespective of the initial guesses.
where J S stands for the saturation current, n is the emission coefficient, V T is the thermal voltage, and J is the diode current. Using the Kirchhoff's second law (V R + V D = V S ) and Ohm's law (V = JR), a root-finding model can be found. The final structure for the model would be as follows: Assuming values of parameters V S , R, n, V T , J S from [38], we obtain the following equation that is nonlinear in the variable J: The approximate solution of the above equation correct to 50 dp is as follows: Table 9 presents numerical simulations for (47) with two different initial conditions of 0.5 and 1.8, under which the smallest absolute error seems to lie under the column of the proposed method PM 10 . Further analysis can easily be conducted by the careful examination of the results tabulated therein. [39].

Example 5. Conversion of Species within a Chemical Reactor
The following nonlinear equation arises in chemical engineering during the conversion of species in a chemical reactor: where x stands for the fractional conversion of the species; thus, it must lie in (0,1). The approximate solution of the above equation correct to 50 dp is as follows: x ≈ 0.75739624625375387945964129792914529342795578042081. Numerical results can be found in Table 10, where it is shown that PM 10 has outperformed all other methods in terms of the absolute errors under consideration under two different initial conditions.  [40]. The two-dimension Bratu system is given by the following partial differential equation: subject to the following boundary conditions The two-dimensional Bratu system has two bifurcated exact solutions for λ < λ c , a unique solution for λ = λ c , and no solutions for λ > λ c . The exact solution to (49) is determined as follows: where θ is an undetermined constant satisfying the boundary conditions and assumed to be the solution of (49). Using the procedure described in [41], one obtains Differentiating (52) with respect to θ and setting dλ dθ = 0, the critical value λ c satisfies By eliminating λ from (52) and (53), we have the value of θ c for the critical λ c satisfying and θ c = 4.798714561. We then obtain λ c = 7.027661438 from (53). Numerical simulations performed in Table 11 show that the proposed three-step method takes a smaller number of iterations and produces considerably smaller absolute errors with a reasonable amount of CPU time. Table 11. Numerical results for problem 6. Now, we consider five different kinds of nonlinear multidimensional equations and numerically solve them with PM 10 , HM 6 , MH 5 , and NR 2 since the methods MH 10 , N M 9 , and WO 8 were either divergent or not applicable on systems of nonlinear equations. For the systems considered, various types of initial guesses are used, and for comparison purposes, the approximate solution and the normed error = ||x n+1 − x n || ∞ having the same tolerance 10 −200 and CPU time are taken into consideration. It can be observed in Tables 12-16 that the smallest possible absolute error is achieved only with the proposed method-that is, PM 10 -in a reasonably affordable time period. Example 7. The nonlinear system of two equations from [3,41] is given as:

Method IG
The exact solution of the system (55) is x = [0, 0] . The numerical results for this system are shown in Table 12 under the proposed PM 10 and other three methods.
where its approximate solution up to 50 dp is shown below: The numerical results for the system (56) are shown in Table 13.

Example 9.
A three-dimensional nonlinear system is taken from [3] as given below: where its approximate solution up to 50 dp is as follows: The numerical results for the system (58) are shown in Table 14. Given below is a nonlinear system of two equations that describe trajectories for the catenary and the ellipse, respectively. We are interested in finding their intersection point that lies in the first quadrant of the cartesian plane. The system has been solved under the proposed PM 10 method and other methods under consideration. The performance of each method is shown in Table 15, whereas an approximate solution up to 50 dp of the system (60) is shown in comparison to the system.
Approximate solution: In this problem, we consider a system developed by Edward Lorenz, who was an American meteorologist studying atmospheric convection around the Earth's surface. Lorenz's nonlinear system is a set of three ordinary differential equations, as given below: In order to study the steady-state behavior of the system (62), we takeẋ 1 (t) =ẋ 2 (t) = x 3 (t) = 0 and a = −1, b = 2, c = 3 to obtain the following nonlinear algebraic system: The approximate solution for the system (63) correct to 50 dp is given as The nonlinear steady-state system (63) has been numerically solved in Table 16.

Concluding Remarks
This research study is based on devising a new, highly effective three-step iterative method with tenth-order convergence. The convergence is proved theoretically via Taylor's series expansion for single and multi-variable nonlinear equations, and the approximate computational order of convergence confirms such findings. Thus, the proposed method PM 10 is applicable not only to single nonlinear equations but also to nonlinear systems. Moreover, dynamical aspects of PM 10 are also explored with basins of attraction that show quite esthetic phase plane diagrams when applied to complex-valued functions, thereby proving the stability of the method when initial guesses are taken within the vicinity of the underlying nonlinear model. Finally, different types of nonlinear equations and systems, including those used in physical and natural sciences, are chosen to be tested with PM 10 and with various well-known optimal and non-optimal methods in the sense of King-Traub. In most of the cases, PM 10 is found to have better results, particularly when it comes to the number of iterations N to achieve required accuracy, ACOC, absolute error, and absolute functional value. It is also worthwhile to note that the proposed method always converges, irrespective of whether the initial guess passed to it lies near to or away from the approximate solution. Hence, PM 10 is a competitive iterative method with tenth-order convergence for solving nonlinear equations and systems.
We understand that methods of very high order are only of academic interest since approximations to the solutions of very high accuracy are not needed in practice. On the other hand, such methods are, to some extent, complicated and do not offer much of an increase in computational efficiency. Moreover, the method proposed in this article lies in the family of methods without memory, requiring the evaluation of three Jacobian matrices, and thereby becomes computationally expensive. To avoid computational complexity, we will propose, in future studies, a modification of the proposed method by replacing the first-order derivative with a suitable finite-difference approximation. In addition, the proposed method will also be analyzed for its semi-local convergence.