The multi-fluid plasma model is derived from moments of the Boltzmann equation for each species. The two-fluid plasma model has electron and ion species. Additional fluids can include neutrals= impurities, or excited ions. Moment equations are truncated, e.g. after five, ten, thirteen, or twenty equations, and are closed with expressions that relate higher moment variables to lower moment variables. Reduction of the two-fluid plasma model to single-fluid MHD is accomplished by applying asymptotic approximations. However, these approximations alter the dispersion relations and can require artificial dissipation for numerical stability. Large mass differences between electrons and ions in the multi-fluid plasma model introduce disparate temporal and spatial scales and require numerical algorithms with sufficient accuracy to capture the multiple scales. Source terms couple the fluids to themselves (interspecies interactions) and to the electromagnetic fields. Interspecies interactions also occur through collisional source terms that account for the direct transfer of momentum and energy. In addition to plasma and electrodynamic physics, the multi-fluid plasma model can capture atomic physics in the form of reaction rate equations for ionization and recombination, which introduce new temporal scales to the plasma dynamics model. The numerical algorithm must treat the inherent stiffness introduced by the multiple physical effects of the model and tightly couple the source terms of the governing equations. The governing equations are expressed in a balance law form and numerical algorithms are developed for their solution. These methods include finite volume and finite element methods with approximate Riemann solvers for the spatial discretization and Runge-Kutta methods for the temporal advance. The algorithms are benchmarked against analytical solutions for wave dispersion and published solutions for the electromagnetic plasma shock problem and the GEM magnetic reconnection problem. In addition, the algorithm has been applied to plasma sheath formation and to study drift turbulence in Z-pinch and field reversed configuration plasmas.