SUMMARY We introduce a new approach to the computation of gravito-elastic free oscillations or normal modes of spherically symmetric bodies based on a spectral element discretization of the radial ordinary differential equations. Our method avoids numerical instabilities often encountered in the classical method of radial integration and root finding of the characteristic function. To this end, the code is built around a sparse matrix formulation of the eigenvalue problem taking advantage of state-of-the-art parallel iterative solvers. We apply the method to toroidal, spheroidal and radial modes and we demonstrate its versatility in the presence of attenuation, fluid layers and gravity (including the purely elastic case, the Cowling approximation, and full gravity). We demonstrate higher-order convergence and verify the software by computing seismograms and comparing these to existing numerical solutions. Finally, to emphasize the general applicability of our code, we show spectra and eigenfunctions of Earth, Mars and Jupiter’s icy moon Europa and discuss the different types of modes that emerge.