The purpose of this work is to develop a general, reliable and efficient algorithm, which is able to deal with multiple reactions in multiphase systems. We selected the method of Lagrange multipliers to minimize the Gibbs energy of the system, under material balance constraints. Lagrange multipliers and phase amounts are the independent variables, whose initialization is performed by solving a subset of the working equations. This initialization is the unconstrained minimization of a convex function and it is bound to converge. The whole solution procedure employs a nested loop with Newton iteration in the inner loop and non-ideality updated in the outer loop, thus giving an overall linear convergence rate. Stability analysis is used to introduce additional phases sequentially so as to obtain the final multiphase solution. The procedure was successfully tested on vapor-liquid equilibrium (VLE) and vapor-liquid-liquid equilibrium (VLLE) of reaction systems. (C) 2017 Elsevier Ltd. All rights reserved.