Systems of nonlinear equations are known as the basis for many models of engineering and data science, and their accurate solutions are very critical in achieving progress in these fields. However, solving a system with multiple nonlinear equations, usually, is not an easy task. Consequently, finding a robust and accurate solution can be a very challenging problem in complex systems. In this work, a novel hybrid method namely Newton-Harris hawks optimization (NHHO) for solving systems of nonlinear equations is proposed. The proposed NHHO combines Newton's method, with a second-order convergence where the correct digits roughly double in every step, and the Harris hawks optimization (HHO) to enhance the search mechanism, avoid local optima, improve convergence speed, and find more accurate solutions. We tested a group of six well-known benchmark systems of nonlinear equations to evaluate the efficiency of NHHO. Further, comparisons between NHHO and other optimization algorithms, including the original HHO algorithm, Particle Swarm Optimization (PSO), Ant Lion Optimizer (ALO), Butterfly Optimization Algorithm (BOA), and Equilibrium Optimization (EO) were performed. The norm of the equation system was calculated as a fitness function to measure the optimization algorithms' performance. A solution with less fitness value is considered a better solution. Furthermore, the experimental results confirmed the superiority of NHHO over the other optimization algorithms, in the comparisons, in different aspects, including best solution, average fitness value, and convergence speed. Accordingly, the proposed NHHO is powerful and more effective in all benchmark problems in solving systems of nonlinear equations compared to the other optimization algorithms. Finally, NHHO overcomes the limitations of Newton's method, including selecting the initial point and divergence problems.INDEX TERMS Hybridization; Harris hawks optimization; Newton's method; optimization; systems of nonlinear equations.