An iterative regularization method is proposed for neutron spectrum unfolding of the liquid organic scintillator EJ-301. The Arnoldi iteration is used to generate Krylov subspaces. GMRES (Generalized Minimal Residual) combined with Tikhonov regularization is used for neutron spectrum updating. The restart method is applied to speed up the convergence as well as non-negativity constraints. The Monte Carlo code Geant4 is applied to simulate the response functions and PHSs (pulse height spectrums) of the neutron spectrums for unfolding. Four typical neutron spectrums are studied, such as the mono-energetic DD fusion neutron spectrum, 241Am-Be neutron spectrum, 252Cf neutron spectrum and Gaussian broadened DT fusion neutron spectrum. Comparison study between the proposed method and the existing methods is also carried out, such as GRAVEL, MLEM, and direct Tikhonov regularization. Unfolding results with acceptable accuracies are obtained with the proposed method, and the proposed method shows better adaptabilities than the other three methods with both the ideal and simulated PHSs in most cases.