In this paper, based on the Ulm-Chebyshev iterative procedure, we present a multi-step Ulm-Chebyshev-like method to solve systems of nonlinear equations $F(x)=0$,

$$\{\begin{array}{l}{\mathbf{y}}_n={\mathbf{x}}_n-B_{n}F({\mathbf{x}}_n),\\ {\mathbf{z}}_n={\mathbf{y}}_n-B_{n}F({\mathbf{y}}_n),\\ {\mathbf{x}}_{n+1}={\mathbf{z}}_n-B_{n}F({\mathbf{z}}_n),\\ {\overline{B}}_n=2B_n-B_n{A}_{n+1}{B}_n,\\ B_{n+1}={\overline{B}}_n+{\overline{B}}_n(2I-A_{n+1}{\overline{B}}_n)(I-A_{n+1}{\overline{B}}_n),n=0,1,2,\dots ,\end{array}$$

where $A_{n+1}$ is an approximation of the derivative $F^{\prime }({\mathbf{x}}_{n+1})$. This method does not contain inverse operators in its expression, and does not require computing Jacobian matrices for solving Jacobian equations. We have proved that the multi-step Ulm-Chebyshev-like method converges locally to the solution with $R$-convergence rate 4 under appropriate conditions. Some applications are given, compared with other existing methods, where the most important features of the method are shown.