We consider two main inverse Sturm–Liouville problems: the problem of recovery of the potential and the boundary conditions from two spectra or from a spectral density function. A simple method for practical solution of such problems is developed, based on the transmutation operator approach, new Neumann series of Bessel functions representations for solutions and the Gelfand–Levitan equation. The method allows one to reduce the inverse Sturm–Liouville problem directly to a system of linear algebraic equations, such that the potential is recovered from the first element of the solution vector. We prove the stability of the method and show its numerical efficiency with several numerical examples.