Recently, a robust new computational technique, based on interval analysis, has been developed for solving the difficult nonlinear problems arising in the modeling of phase behavior. This technique can be used, with mathematical and computational guarantees of certainty, to find the global optimum of a nonlinear function or to enclose any and all roots of a system of nonlinear equations. As shown in the applications here to phase stability analysis and to the location of homogeneous azeotropes, it provides a method that can guarantee that the correct result is found; thus, eliminating computational problems that may potentially be encountered with currently available techniques. The method is initialization independent; it is also model independent, straightforward to use, and can be applied in connection with any equation of state or activity coefficient model.