The multi-fluid plasma model only assumes local thermodynamic equilibrium within each fluid, e.g. ion and electron fluids for the two-fluid plasma model. Derivation of the MHD model involves several asymptotic and simplifying assumptions that can limit its applicability. Therefore, the two-fluid plasma model more accurately represents the appropriate physical processes. Physical parameters indicate the importance of the two-fluid effects: electron to ion mass ratio m e /m i , ion skin depth δ i , and ion Larmor radius r L . The MHD model assumes m e /m i = 0, δ i = 0, and r L = 0. Asymptotic approximations of the two-fluid model, Hall-MHD, has an unbounded Whistler wave that requires artificial dissipation. No unbounded waves exist in the two-fluid model. An algorithm is developed for the simulation of plasma dynamics using the two-fluid and multi-fluid plasma models. The algorithm implements a discontinuous Galerkin method that uses an approximate Riemann solver[1] to compute the fluxes of the fluids and electromagnetic fields at the computational cell interfaces. The two-fluid plasma model has time scales on the order of the electron and ion cyclotron frequencies, the electron and ion plasma frequencies, the electron and ion sound speeds, and the speed of light. Such disparate time scales motivate a semi-implicit time-stepping scheme to overcome the severe time step restrictions of explicit schemes. The algorithm is validated with several test problems including the GEM challenge magnetic reconnection problem [2] and the generation of dispersive plasma waves which are compared to analytical dispersion diagrams. The algorithm is applicable to study advanced physics calculations of plasma dynamics including magnetic plasma confinement and astrophysical plasmas. Three-dimensional solutions of the Z-pinch and the field reversed configuration (FRC) magnetic plasma confinement configurations are presented.