The multi-fluid plasma model is derived from moments of the Boltzmann equation and typically has two fluids representing electron and ion species. Large mass differences between electrons and ions introduce disparate temporal and spatial scales and require a numerical algorithm with sufficient accuracy to capture the multiple scales. Source terms of the multi-fluid plasma model couple the fluids to themselves (interspecies interactions) and to the electromagnetic fields. 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. A discontinuous Galerkin method is implemented for the spatial representation. Time integration is investigated using explicit, implicit, semi-implicit methods. Semi-implicit treatment is accomplished using a physics-based splitting. The algorithm is applied to study drift turbulence in field reversed configuration plasmas to illustrate the physical accuracy of the model. The algorithm is also applied to plasma sheath formation which demonstrates Langmuir wave propagation.