The stable computation of linearly constrained, multiphase, chemical equilibrium compositions is an important topic for a wide range of industrial and academic applications. Numerous computational methods have been developed to solve such problems which are, in general, susceptible to failure under certain conditions due to numerical stiffness. In this work, we present a Gibbs function continuation method for linearly constrained multiphase equilibrium calculations. The method converts the nonlinear Element Potential Equations - derived from the minimization of the mixture Gibbs free energy using the Lagrange multiplier technique - into an initial value problem which can be stably integrated through the use of a property of linear least squares solutions. The stability and convergence properties of the proposed method are derived and it is shown that the single phase method arises as a special case of the multiphase algorithm when only one phase is considered. Two test cases are presented to clarify and demonstrate the accuracy and robustness of the method. (C) 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.