We develop a lattice Boltzmann equation method for simulating multi-phase immiscible fluid flows with variation of density and viscosity, based on the model proposed by Gunstensen {\em et al} for two-component immiscible fluids. The numerical measurements of surface tension and viscosity agree well with theoretical predictions. Several basic numerical tests, including spinodal decomposition, two-phase fluid flows in two-dimensional channels and two-phase viscous fingering, are shown in agreement of experiments and analytical solutions.