The paper deals with the numerical solution of inverse Sturm-Liouville problems with unknown potential symmetric over the interval [0, ]. The proposed method is based on the use of a family of Boundary Value Methods, obtained as a generalization of the Numerov scheme, aimed to the computation of an approximation of the potential belonging to a suitable function space of finite dimension. The accuracy and stability properties of the resulting procedure for particular choices of such function space are investigated. The reported numerical experiments put into evidence the competitiveness of the new method.