  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.
