Finding all optimal solutions for a metabolic model is the challenge of metabolic modeling, but there is no practical algorithm for large scale models. A two-phase algorithm is proposed here to systematically identify all optimal solutions. In phase 1, the model is reduced using the FVA approach; in phase 2, all optimal solutions are searched by the addition of a binary variable to convert the model to an MILP problem. The proposed approach proved itself to be a more tractable method for large scale metabolic models when compared with the previously proposed algorithm. The algorithm was implemented on a metabolic model of Escherichia coli (iJR904) to find all optimal flux distributions. The model was reduced from 1076 to 80 fluxes and from 998 to 54 equations and the MILP problem was solved, resulting in 30,744 various flux distributions. For the first time, this number of optimal solutions has been reported. (C) 2014 Elsevier Ltd. All rights reserved.