At fixed $l$, for each $n$ above the label of first electronic level (ground state), the corresponding eigenstate must be orthogonal to the previous ones. It is this constraint which forces the wavefunction to give larger weight on distances in average larger than the previous $n$ states, thus justifying a larger energy.
More formally, looking at the solution of a Schrödinger-like equation as a variational problem, the ground state is obtained by minimizing
$$
\left<\psi_0 \left| H \right| \psi_0 \right > - \lambda (\left<\psi_0 \left|\right. \psi_0 \right > -1 )
$$
with respect to the state $\left|\psi_0 \right >$, $\lambda$ being the Lagrange multiplier which takes care of the normalization constraint and which corresponds to the energy. The extremum condition with respect to $\left|\psi_0 \right >$ provides the time independent Schrödinger-like equation
$$
\left. \left. H \right| \psi_0 \right >= \lambda \left|\psi_0 \right >
$$
I am using the term "Schrödinger-like" because the same considerations apply to a pure Schrödinger equation but also to non-linear equations corresponding to all kind of self-consistent field approsimations, including the Kohn-Sham equations coming from Density Functional Theory for the electronic states.
The next-above-the-minimum-energy state will be the solution of a variational problem for a new state $\left|\psi_1 \right >$, where an additional constraint is added, thus an additional term proportional to the scalar product $\left<\psi_1 \left|\right. \psi_0 \right >$ appears. The resulting Schrödinger-like equation will contain a term proportional to $\left|\psi_0 \right >$, i.e. the ground state wavefunction acts as an additional interaction term which, in coordinate representation, forces the state $\left|\psi_1 \right >$ to have a nodal surface and to assign a higher weight to larger distances. By construction, this state must have an energy larger than the ground state. The analysis of the effect of the extra term allows to rationalize the origin of such larger energy.