In this thesis, we develop a general framework for the lattice Boltzmann method to simulate multiphase systems with an arbitrary number of components. Theoretical expectations are easily visualized for binary mixtures, so we focus on characterizing the performance of the method by numerically minimizing the free energy of a binary van der Waals mixture to generate phase diagrams. Our phase diagrams contain very intriguing features that are not well-known in today’s physics community but were understood by van der Waals and his colleagues at the turn of the 20th century. Phase diagrams and lattice Boltzmann simulation results are presented in a density-density plane, which best matches with LB simulations performed at constant volume and temperature. We also demonstrate that the algorithm provides thermodynamically consistent results for mixtures with larger numbers of components and high density ratios. All of the theoretical phase diagrams are recovered well by our lattice Boltzmann method.