A new algorithm to determine liquid-liquid phase equilibrium at the global minimum of the Gibbs free energy is presented, the optimization results are then used by Newton-Raphson method to obtain the final accurate phase equilibrium compositions. The tangent plane distance function criterion is used as the sufficient and necessary condition for the global stability analysis, its global minimum is obtained by the simulated annealing algorithm in order to overcome the multiple local minimum problem. The Gibbs free energy functions for liquid-liquid equilibrium systems are formulated by NRTL or UNIQUAC activity coefficient equation. Aiming to overcome the nonconvex problem of the Gibbs free energy function the simulated annealing algorithm is used for achieving a near true equilibrium solution, this solution then is used as the initial values of the Newton-Raphson method. Our algorithm was implemented and tested for six examples with three components and up to three liquid phases. In all of the examples, results for liquid-liquid equilibrium calculations and global stability analysis were found to be reliable. The main advantage of our algorithm is its simplicity and that it is model independent.