We develop a numerical algorithm for the solution of the Sturm-Liouville differential equation governing the stationary radial oscillations of nonrotating compact stars. Our method is based on the Numerov's method that turns the Sturm-Liouville differential equation in an eigenvalue problem. In our development we provide a strategy to correctly deal with the star boundaries and the interfaces between layers with different mechanical properties. Assuming that the fluctuations obey the same equation of state of the background, we analyze various different stellar models and we precisely determine hundreds of eigenfrequencies and of eigenmodes. If the equation of state does not present an interface discontinuity, the fundamental radial eigenmode becomes unstable exactly at the critical central energy density corresponding to the largest gravitational mass. However, in the presence of an interface discontinuity, there exist stable configurations with a central density exceeding the critical one and with a smaller gravitational mass.