The objective of this study is to explore the application of high-order numerical methods in ionosphere-plasmasphere modeling. Specifically, the nodal discontinuous Galerkin (DG) method is chosen to solve the multifluid dynamical equations along the magnetic field lines. A general curvilinear magnetic field-line-following coordinate system is also used in the model. Numerical simulations with different combinations of number of elements (K) and polynomial orders (N) show the converging results, indicating the robustness of the algorithms and implementation. The model also captures the dawn terminator effect very well in the He⁺ field.