Let us rewrite Schrödinger's equation as
\[ \Psi'' +\frac{2m}{\hbar^2}(E-U) \Psi = 0, \quad k^2 = \frac{2m}{\hbar^2}(E-U)\]Its solutions are the sum of harmonic functions
\[ \Psi = A \cos kx + B \sin kx.\]Taking into account the boundary conditions $\Psi(0)=\Psi(a)=0$, we have a non-trivial solution when $A=0$ and $ka = \pi n$. Taking into account $U=0$, we have
\[\Psi_n = B \sin \frac{\pi n x}{a}, \quad E_n = \frac{\hbar^2 \pi^2 n^2}{2m a^2}\]
The first electron fills the level with energy $E_1$, the second with energy $E_2$, and so on. Therefore,
\[E = \frac{\hbar^2 \pi^2}{2ma^2} \sum_1^Nn^2 =\frac{\hbar^2 \pi^2}{2ma^2} \frac{N(N+1)(2N+1)}{6} \simeq \frac{\hbar^2 \pi^2}{2ma^2} \frac{N^3}{3} \]It is correct to take into account the spin of the electron in this calculation. Then the electrons fill the levels not one by one, but in pairs. And the answer then (in the approximation $N \gg 1$)
\[E = \frac{\hbar^2 \pi^2}{2ma^2} \frac{N^3}{24}\]
Let us find the force $F$ using the method of virtual displacements. When the width of the well changes by $da$, the work of forces $F$ (on electrons) and the change in energy must be compensated:
\[2F \, da - \frac{\hbar^2\pi^2 N^3}{3ma^3} \, da = 0 \quad \Rightarrow \quad F = \frac{\hbar^2 \pi^2 N^3}{6ma^3}\]