Inertial waves, which are predominantly driven by the Coriolis force, likely play an important role in solar dynamics, and, additionally, they provide a window into the solar subsurface. The latter allows us to infer properties that are inaccessible to the traditional technique of acoustic wave helioseismology. Thus, a full characterization of these normal modes holds the promise of enabling investigations into solar subsurface dynamics. In this work, we develop a spectral eigenvalue solver to model the spectrum of inertial waves in the Sun. We model the solar convection zone as an anelastic medium, and solve for the normal modes of the momentum and energy equations. We demonstrate that the solver can well reproduce the observed mode frequencies and line widths, not only of sectoral Rossby modes, but also of recently observed high-frequency inertial modes. In addition, we believe that the spectral solver is a useful contribution to the numerical methods of modeling inertial modes on the Sun.