Several engineering applications need a robust method to find all the roots of a set of nonlinear equations automatically. The proposed method guarantees monotonous convergence, and it can determine whole submanifolds of the roots if the number of unknowns is larger than the number of equations. The critical steps of the multidimensional bisection method are described and possible solutions are proposed. An efficient computational scheme is introduced. The efficiency of the method is characterized by the box-counting fractal dimension of the evaluated points. The multidimensional bisection method is much more efficient than the brute force method. The proposed method can also be used to determine the fractal dimension of the submanifold of ...