## 1 Introduction  
An often discussed problem in many applications of science and technology is to find a real zero of a system of nonlinear equations \(F(x)=0\), where \(F(x)=(f_{1}(x),f_{2}(x),...,f_{n}(x))^{T}\), \(x=(x_{1},x_{2},...,x_{n})^{T}\)\(f_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R},\forall i=1,2,\ldots,n\) and \(F:D\subset\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\) is a smooth map and \(D\) is an open and convex set, where we assume that \(\alpha=(\alpha_{1},\alpha_{2},...,\alpha_{n})^{T}\) is a zero of the system and \(x^{(0)}=\left(x_{1}^{(0)},x_{2}^{(0)},...,x_{n}^{(0)}\right)^{T}\) is an initial guess sufficiently close to \(\alpha\). For example, problems of the above type arise while solving boundary value problems for differential equations. The differential equations are reduced to system of nonlinear equations, which are in turn solved by the familiar [PERSON]'s iteration method having convergence order two [1]. The [PERSON] method (\(2^{nd}\)_NM_) is given by  
\[x^{(k+1)}=G_{2^{nd}NM}(x^{(k)})=x^{(k)}-u(x^{(k)}),\quad u(x^{(k)})=[F^{\prime} (x^{(k)})]^{-1}F(x^{(k)}) \tag{1}\]  
[PERSON] [2] has proposed a third order iterative method called Harmonic Mean [PERSON]'s method for solving a single nonlinear equation. Analogous to this method [2], we consider the following extension to solve a system of nonlinear equation \(F(x)=0\), henceforth called as \(3^{rd}HM\):  
\[x^{(k+1)}=G_{3^{rd}HM}(x^{(k)})=x^{(k)}-\frac{1}{2}\left([F^{\prime}(x^{(k)})] ^{-1}+[F^{\prime}(x^{(k)}-u(x^{(k)}))]^{-1}\right)F(x^{(k)}) \tag{2}\]  
We note that \(\frac{1}{2}\left([F^{\prime}(x^{(k)})]^{-1}+[F^{\prime}(x^{(k)}-u(x^{(k)}))]^{ -1}\right)\) is the average of the inverses of two Jacobians. In general, such third order methods free of second derivatives like Equation (2) can be used for solving system of nonlinear equations. These methods require one function evaluation and two first order Frechet derivative evaluations. The convergence analysis of a few such methods using point of attraction theory can be found in [3]. This \(3^{rd}HM\) method is more efficient than [PERSON]'s method because it does not require the evaluation of a third order tensor of \(n^{3}\) values while finding the number of function evaluations. Furthermore, the \(3^{rd}HM\) methods are less efficient than two-step fourth order [PERSON]'s method (\(4^{th}NR\))  
\[x^{(k+1)}=G_{4^{th}NR}(x^{(k)})=G_{2^{nd}NM}(x^{(k)})-F^{\prime}(G_{2^{nd}NM}( x^{(k)}))^{-1}F(G_{2^{nd}NM}(x^{(k)})) \tag{3}\]  
which was recently rediscovered by [PERSON] _et al._[4] using the variational iteration technique. Recently [PERSON] _et al._[5] developed the fourth order method, which is given by  
\[ \begin{split}& x^{(k+1)}=G_{4^{th}SGS}(x^{(k)})=x^{(k)}-W(x^{(k)})u (x^{(k)}),\\ & W(x^{(k)})=-\frac{1}{2}I+\frac{9}{8}[F^{\prime}(y(x^{(k)}))]^{- 1}F^{\prime}(x^{(k)})+\frac{3}{8}[F^{\prime}(x^{(k)})]^{-1}F^{\prime}(y(x^{(k )}))\\ & y(x^{(k)})=x^{(k)}-\frac{2}{3}u(x^{(k)})\end{split}  \tag{4}\]  
[PERSON] _et al._[6] presented a sixth order method, which is given by  
\[ \begin{split}& x^{(k+1)}=G_{6^{th}CHMT}(x^{(k)})=z(x^{(k)})-[F^{ \prime}(x^{(k)}-u(x^{(k)}))]^{-1}F(z(x^{(k)}))\\ & z(x^{(k)})=x^{(k)}-u(x^{(k)}-\Big{[}2I-F^{\prime}(x^{(k)})^{-1} F^{\prime}(x^{(k)}-u(x^{(k)}))\Big{]}[F^{\prime}(x^{(k)})]^{-1}F(x^{(k)}-u(x^{(k)})) \end{split}  \tag{5}\]  
Recently, an improved fourth order version from a third order method for solving a single nonlinear equation is found in [7]. In the current paper, similar to the method found in [7], a multivariate version having fourth order convergence is proposed. The rest of this paper is organized as follows. In Section 2, we present a new algorithm (optimal) that has fourth order convergence by using only three function evaluations and a multi-step version with order \(2r+4\), where \(r\) is a positive integer and \(r\geq 1\) for solving systems of nonlinear equations. In Section 3, we study the convergence analysis of the new methods using the point of attraction theory.