In this paper, we present a neutron spectrum unfolding technique based on a modification of the least-squares method. The main innovation is the use of a Krylov subspace iteration method to solve the least-squares normal equations. This method was employed because it performs better on ill-conditioned systems of linear equations as compared with standard direct-solution methods. Three different least-squares solution techniques are compared and evaluated in terms of (i) accuracy in the prediction of the energy spectrum, (ii) computational efficiency, and (iii) robustness to noise. The unfolding is performed on measured pulse-height distributions as well as pulse height distributions generated with the Monte Carlo code MCNP-PoliMi. Using this code, neutron energy depositions on the constituents of the scintillator are individually tracked, and the light output generated at each interaction is suitably accounted for. This procedure allows for a very accurate simulation of the liquid scintillator detector response. The precise knowledge of the neutron energy spectrum provides information not only about the presence or absence of fissile material, but also about the characteristics of the material. We show that the proposed technique performs well in the unfolding of neutron pulse-height distributions from Monte Carlo simulations, and fairly well for a measured distribution from a Cf-252 neutron source.