I don't know how to solve this problem just like that, but I'd suggest the following approach. Because we know that the two electrons are confined to a finite interval, we know that their respective hilbert spaces in the space representation take on the following form:
psi_n(x) = N sin(n pi x / L)
The product hilbert space then has a basis psi_{n,m}(x,y) = psi_n(x) cross psi_m(y)
which we can represent by the dirac ket |n,m>.
If we ignore the fermionic part, in that basis, we can work out the hamiltonian from the coulomb interaction:
E_{k,l} + E_{m,n} + <k,l | 1/|x-y| | m,n>
the last term is an integral which I think can be worked out, and the first two E-terms are the free energy terms for the wave functions psi_n etc...
Now the final problem - which I don't know how to tackle - is how to diagonalise this infinite matrix.
If we want to introduce the fermionic character of the electrons, the modification is essentially that we should only consider anti-symmetric combinations (a redefinition of |k,l> as a function of the sine functions).
cheers,
Patrick.