Lattice Boltzmann methods (LBM) have several features that make them attractive for computation of large fluid mechanics problems in complex geometries. For example, in comparison with the conventional projection methods that are widely used in computational fluid dynamics, LBM does not require solution of a pressure Poisson equation which usually accounts for similar to 80% of the computational time required per time step, and it is also local, making parallelization straightforward, and run time on distributed memory architectures linearly scalable. LBM has therefore been considered for direct simulation of multiphase flows, but its application has been limited because of instabililties that arise when the viscosities become small and/or the density mismatch between the fluids is large. Current approaches to computation of multiphase systems using LBM are reviewed, and new approaches based on multiple relaxation times and a regularization procedure which maintain stability at low viscosities are discussed and a technique using time-splitting, that alleviates the density ratio constraint, is proposed. Applications of the LBM to magnetohydrodynamic (MHD) multiphase flows will be discussed.