Why does this spatially symmetric wavefunction have lower energy than its anti-symmetric counterpart?

Here's a hand-wavy argument that runs like this. In the Schrodinger equation: $$ -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x)+V(x)\psi(x)=E\psi(x) $$ the second derivative term $-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x)$ plays the role of the kinetic energy. This is basically the curvature of the wavefunction.

For, for fixed potential, smaller curvature, i.e. smaller values of $-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x)$, will produce a smaller energy $E$.

Next, note that odd solutions are necessarily $0$ at the origin, and that the solutions must all to go zero as $x\to \pm \infty$.

Given this, an even (or symmetric) solution can have very low curvature since it can "spread out" over the well without having to go through $0$ at the origin.

Indeed, a property of the solutions to the Schrodinger equation is that the energy of the solution must increase with the number of nodes. Following Albert Messiah's book on Quantum Mechanics, the oscillatory theorem states that, if \begin{equation} \hat H\psi _{1}(\xi )=E_{1}\psi _{1}(\xi )\, ,\qquad \hat H\psi_{2}(\xi)=E_{2}\psi _{2}(\xi)\, ,\qquad E_{1}<E_{2}, \end{equation} then $\psi _{2}(\xi )$ has more zeroes (nodes) than $\psi _{1}(\xi ).$ To show this, manipulate the Schrodinger equation so that \begin{equation} \frac{1}{\psi _{1}(\xi )}\frac{\hbar ^{2}}{2m}\frac{d^{2}}{d\xi ^{2}}\psi _{1}(\xi )+E_{1}=V(\xi )=\frac{1}{\psi _{2}(\xi )}\frac{\hbar ^{2}}{2m}\frac{% d^{2}}{d\xi ^{2}}\psi _{2}(\xi )+E_{2}, \end{equation} which one then reorganizes as \begin{eqnarray} \psi _{2}(\xi )\frac{d^{2}}{d\xi ^{2}}\psi _{1}(\xi )-\psi _{1}(\xi )\frac{ d^{2}}{d\xi ^{2}}\psi _{2}(\xi ) &=&\frac{d}{d\xi }\left( \psi _{2}(\xi ) \frac{d}{d\xi }\psi _{1}(\xi )-\psi _{1}(\xi )\frac{d}{d\xi }\psi _{2}(\xi )\right) \nonumber \\ &=&\frac{\hbar ^{2}}{2m}\left( E_{2}-E_{1}\right) \psi _{1}(\xi )\psi _{2}(\xi ). \end{eqnarray} Integrating from $\xi =a$ to $\xi =b,$ we find \begin{eqnarray} &&\frac{\hbar ^{2}}{2m}\left( E_{2}-E_{1}\right) \int_{a}^{b}d\xi \psi _{1}(\xi )\psi _{2}(\xi ) \nonumber \\ &&\qquad = \left( \psi _{2}(b)\frac{d}{d\xi }\psi _{1}(b)-\psi _{1}(b)\frac{d}{d\xi }\psi _{2}(b)\right)\nonumber \\ &&\qquad\qquad -\left( \psi _{2}(a)\frac{d}{d\xi }\psi _{1}(a)-\psi _{1}(a)\frac{d}{d\xi } \psi _{2}(a)\right) \end{eqnarray} If one chooses the points $a,$ $b$ to be consecutive zeroes of $\psi _{1},$ and if $\psi _{1}$ is continuous, then $\psi _{1}$ does not change sign over the interval $\left[ a,b\right] .$ We can assume, without loss of generality, that $\psi _{1}$ is positive on this interval, so that $\psi _{1}^{\prime }(a)\geq 0$ and $\psi _{1}^{\prime }(b)\leq 0.$ Then, using these with $\psi _{1}(a)=\psi _{1}(b)=0,$ we obtain \begin{equation} \frac{\hbar ^{2}}{2m}\left( E_{2}-E_{1}\right) \int_{a}^{b}d\xi \psi _{1}(\xi )\psi _{2}(\xi )=\psi _{2}(b)\frac{d}{d\xi }\psi _{1}(b)-\psi _{2}(a)\frac{d}{d\xi }\psi _{1}(a). \end{equation}

Assume now that $\psi _{2}(\xi )$ is always positive on $\left[ a,b\right] .$ Then the integrand is positive, so the integral will be positive (the area under a positive curve is necessarily positive!) On the right hand side, however, the first term is negative by virtue of $\psi _{1}^{\prime }(b)\leq 0$ while the second is positive. Overall, the right hand side is therefore negative, while the left side is negative, a contradiction : $\psi _{2}(\xi )$ cannot be always positive on $\left[ a,b\right] ,$ meaning that between two consecutives zeroes of $\psi _{1}(\xi )$, the function $\psi _{2}(\xi )$ must itself go through zero. In other words, $\psi _{2}(\xi )$ has more zeroes than $\psi _{1}(\xi ).$ This proves the theorem.

This result, applied to your problem, shows that the most symmetric solution, which has no nodes, must have the lowest energy.

The theorem is also beautifully illustrated with the harmonic oscillator wave function. By inspection, the wave function for the lowest energy has no zeroes, the wave functions for the next lowest energies have one, two, three, etc. nodes, the number of nodes increasing with the energy. It is also easily verified that, between two nodes of a solution, another solution of higher energy will go through zero at least once.

BTW this also works for the Coulomb potential provided $V(r)$ is the effective potential, which includes the centrifugal term. For given fixed $\ell$ the number of nodes also increases with energy, as expected.