We introduce a methodology to solve nonlinear systems of equations with bound constraints, and two or more roots. In order to find more than one root, this methodology uses an appropriate global optimization algorithm together with a polarization technique. Polarization modifies (successively after the determination of a new root) a merit function associated to the nonlinear system, creating repulsive regions in the neighborhood of the previous roots. We applied successfully the proposed methodology to solve hard nonlinear systems, including the double retrograde vaporization problem for binary systems of methane + n-butane. (C) 2009 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.