A method is described for solving the energy-dependent neutron diffusion equation by first factorizing the flux into a spatial shape function with weak energy dependence and a spectral function, then developing coupled equations for these two functions which must be solved iteratively. Numerical procedures used to solve these equations combine internally, and in a self-consistent fashion, a fine-group spectrum calculation with a broad-group spatial calculation. Numerical examples, based on representative fast-reactor models, are presented to demonstrate that this space-energy factorization method constitutes an accurate and economical approximation.