If we consider a cubic shaped quantum box and put a charged particle in the ground state there, we can compute the electrostatic potential at any point inside the cube using the integral:
$$I=\int_{-1}^1 \int_{-1}^1 \int_{-1}^1 \frac{\cos^2 (\frac{\pi}{2} x) \cos^2 (\frac{\pi}{2} y) \cos^2 (\frac{\pi}{2} z) ~dx dy dz}{\sqrt{(x-X)^2+(y-Y)^2+(z-Z)^2}}$$
Then $\phi(aX,aY,aZ)= -\frac{q}{4 \pi \varepsilon_0 \varepsilon a} I$ where $a$ is the cube side length and $q$ the charge.
It's a very simple problem, however I wanted to know if it's possible to evaluate this kind of integral analytically.
I don't believe we can do it by "brute force", i.e. integrating w.r.t. each coordinate in turn. Which is why I had a thought to use the divergence theorem. Then we could reduce the triple integral to 6 double integrals over the faces of the cube.
But I don't know how to represent the integrated function as the divergence of some field. I would appreciate some ideas or references.