An approximation that seems to work well for the multi-electron case is the Hartree-Fock method.
In Hartree-Fock, we assume the mean-field approximation. Each electron feels the repulsion from other electrons based on their average, not instantaneous, positions. (This assumption prevents Hartree-Fock from predicting van der Waals forces.)
We thus modify the hydrogen Hamiltonian by introducing two new operators. One is the average Coulombic repulsion between electrons, and the other is the exchange interaction. However, because we're using the average position of the electrons, then for our spherical atom these operators don't have an angular dependence. Thus the spherical harmonics are still separable as in the hydrogen case, so roughly the shape of the orbitals must remain the same. The only part that can change is the radial part of the wavefunction. Doing the calculations, you'll see that the radial part of the wavefunctions are squeezed or stretched a little bit due to Coulombic repulsion and the exchange interaction between electrons, and the increased Coulombic attraction to the nucleus. But as Wikipedia says, qualitatively they don't change much until you introduce multiple atoms.
Without the mean-field approximation, I suppose even the angular shape would change, but that's beyond me.