Electric potential of a spheroidal gaussian

I'm looking for results that compute the electrostatic potential due to a spheroidal gaussian distribution. Specifically, I'm looking for solutions of equations of the form $$\nabla^2\Phi=N\exp\left({-\rho^2/2\sigma_r^2}\right)\exp\left({-z^2/2\sigma_\ell^2}\right),$$ where $\rho^2=x^2+y^2$, possibly with some polynomial factor of the form $\rho^{|m|}\cos(m\phi)z^n$ in front.

This has an easy solution when the two variances $\sigma_r$ and $\sigma_\ell$ are equal, in which case there is spherical symmetry applies and Gauss's law readily yield $\Phi$ in terms of error functions.

However, I'm having trouble finding results for the general case. Using spheroidal coordinates doesn't seem to help, as the constant-$\mu$ surfaces are confocal spheroids that taper out to spheres as they get larger, instead of maintaining their ellipticity as the charge density does.

This looks standard enough that it ought to have been done before (right?), but it is also messy enough that I'm unsure it made it to any textbook. Has anyone seen the likes of this before?

-
Just took a quick look, but why doesn't separation of variables work? –  BebopButUnsteady Mar 18 '13 at 23:40
Using cylindrical harmonics the $z$ integral comes out to $\int_{-\infty}^\infty \exp\left(-\frac{z'^2}{2\sigma_\ell}-k|z-z'|\right)\mathrm{d}z'$ which is doable, in terms of error functions, but then chokes the $k$ integral. –  Emilio Pisanty Mar 19 '13 at 0:05