What you probably want is the Euler-Chebychev method for functions $F:\Bbb R^n\to\Bbb R^n$, a multi-dimensional generalization of the Halley method.
You know that $s=-F'(x)^{-1}F(x)$ gives $F(x+s)=O(s^2)$ and thus $(x+s)-x^*=O(s^2)$. Now add a further refinement $v=O(s^2)$ based on the values at the point $x$ and disregard any terms $O(s^3)$,
\begin{align}
F(x+s+v)&=F(x)+F'(x)(s+v)+\frac12F''(x)[s+v,s+v]+O((s+v)^3)\\
&=F'(x)v+\frac12F''(x)[s,s]+O(s^3)
\end{align}
Thus with $v=-\frac12F'(x)^{-1}F''(x)[s,s]$ (which is indeed $O(s^2)$) you obtain a third order method,
\begin{align}
s_n&=-F'(x_n)^{-1}F(x_n)\\
x_{n+1}&=x_n+s_n-\frac12F'(x_n)^{-1}F''(x_n)[s_n,s_n]
\end{align}
For non-square systems replace the inverse matrix by the Moore-Penrose pseudo-inverse of $F'(x_n)$ to get cubic convergence to a point of minimal error.