Calculation of phase equilibria is a crucial part in modeling and simulation of many chemical engineering processes. This is in particular true for systems that exhibit two liquid phases, since methods based on local solvers (e.g. Newton type) may converge to trivial solutions. In this contribution, a new method to test for phase stability and to calculate the globally stable phases for systems that exhibit two liquid phases is presented. The method is based on a preprocessing step, in which all heterogeneous regions are determined. The actual phase stability test and flash calculation is performed using a homotopy continuation algorithm. The method is applied to the LLE flash, the bubble and dew point problem, and to a general (VLLE) p,T-flash. The algorithm is shown to be very reliable while its computational demand is small compared with other methods. Hence, the method seems to qualify for applications in simulation and optimization where flash calculations are performed many times.