We present a variational algorithm for solving the classical inverse Sturm-Liouville problem in one dimension when two spectra are given. All critical points of the least-squares functional are at global minima, which suggest minimization by a (conjugate) gradient descent algorithm. Numerical examples show that the resulting algorithm works quite reliably without tuning for particular cases, even in the presence of noise. With the right choice of parameters, the minimization procedure converges very quickly.