Density functional theory (DFT) is the de facto method for the electronic structure of weakly correlated systems. But for strongly correlated materials, common density functional approximations break down. Here, we derive a many-electron expansion (MEE) in DFT that accounts for successive one-, two-, three-, ... particle interactions within the system. To compute the correction terms, the density is first decomposed into a sum of localized, nodeless one-electron densities (rho(i)). These one-electron densities are used to construct relevant two- (rho(i) + rho(j)), three- (rho(i) + rho(j) +rho(k)), ... electron densities. Numerically exact results for these few-particle densities can then be used to correct an approximate density functional via any of several many-body expansions. We show that the resulting hierarchy gives accurate results for several important model systems: the Hubbard and Peierls-Hubbard models in 1D and the pure Hubbard model in 2D. We further show that the method is numerically convergent for strongly correlated systems: applying successively higher order corrections leads to systematic improvement of the results. MEE thus provides a hierarchy of density functional approximations that applies to both weakly and strongly correlated systems.