Recently, we have seen good progress in our capability to simulate complex electromagnetic systems. However, still there exist many challenges that have to be tackled in order to push limits restricting the field of computational electromagnetics upward. One of these challenges is the limitations in the available computational resources. Over several decades, the traditional computational methods, such as finite difference, finite element, and finite volume methods, have been extensively applied in the field of electromagnetics. On the other hand, the spectral element method (SEM) has been recently utilized in some branches of electromagnetics as waveguides and photonic structures for the sake of accuracy. In this paper, the numerical approximation to the set of the partial differential equations governing a typical magnetostatic problem is presented by using SEM for the first time to the best of our knowledge. Legendre polynomials and Gauss-Legendre-Lobatto grids are employed in the current study as test functions and meshing of the elements, respectively. We also simulate a magnetostatic problem in order to verify the SEM formulation adapted in the current study.