This could also be solved numerically. If the centers of the boxes are at, say, points ##c_1## and ##c_2## on the ##x##-axis, and the length of both boxes is ##L##, the potential ##V(x_1 ,x_2 )## could be approximated with a function that jumps to a very large but finite value when one of the electrons crosses the wall of a box and where the Coulomb term is proportional to ##\frac{1}{|x_1 - x_2 |+\epsilon}## with ##\epsilon## a very small positive number that prevents the potential from being infinite in any situation. Then we could choose an antisymmetrized initial state like
##\Psi (x_1 ,x_2 ,t_0 ) = C\left[\exp\left(-10\left(\frac{x_1 - c_1}{L}\right)^2\right)\exp\left(-10\left(\frac{x_2 - c_2}{L}\right)^2\right) - \exp\left(-10\left(\frac{x_2 - c_1}{L}\right)^2\right)\exp\left(-10\left(\frac{x_1 - c_2}{L}\right)^2\right)\right]##,
and propagate it forward with an imaginary time coordinate as in these blog posts of mine:
https://physicscomputingblog.com/20...art-5-schrodinger-equation-in-imaginary-time/
https://physicscomputingblog.com/20...ution-of-pdes-part-7-2d-schrodinger-equation/ .
Then the state would decay to an approximate ground state of the potential of the original problem (however, if this is done for a long enough time interval, the electrons would tunnel between boxes and end up in an "actual" ground state where there would be possible to find two electrons in the same box with a nonzero probability).