$V_0=\frac{\rho_0}{4\pi \epsilon_o}\iiint_0^{\infty}\frac{e^{-(x^2+2y^2+2z^2)}}{\sqrt{x^2+y^2+z^2}} dxdydz$ Here is what i'm given:
"The so called Schwinger parametrization is based on identities like the following:
$$\frac{1}{\sqrt{R}}=\frac{2}{\sqrt{\pi}}\int_0^{\infty}e^{-Ru^2}du$$
The first octant of a three-dimensional coordinate system is filled with continuously distributed charge of density $\rho(x,y,z)=\rho_0 e^{-(x^2+2y^2+2z^2)}$. Use the above relation to evaluate the electrostatic potential in the origin, which is given by the integral:
$$V_0=\frac{\rho_0}{4\pi \epsilon_o}\iiint_0^{\infty}\frac{e^{-(x^2+2y^2+2z^2)}}{\sqrt{x^2+y^2+z^2}} dxdydz$$
Hint: The Schwinger parametrization converts the triple integral into a quadrupple integral. Three of the four integrals are Gaussian and thus can be easily evaluated. Use a computer algebra system to evaluate the last one-dimensional integral that you obtain."
This is what i've tried (but i am probably completely off base):
This problem looks like an obvious candidate for spherical coordinates so using the the identities $x=r\sin\phi\cos\theta$, $y=r\sin\phi\sin\theta$, $z=r\cos\phi$ i have tried to simplify the terms in the exponent as such:
$$Exp[-(x^2+2y^2+2z^2)]=Exp[-r^2(\sin^2\phi\cos^2\theta + 2\sin^2\phi\sin^2\theta + 2\cos^2\phi)]=Exp[-r^2(\frac{1-\cos2\phi}{2}\frac{1+\cos 2\theta}{2})+(1-\cos 2\phi)(1-\cos 2\theta)+(1+\cos 2\phi))]=Exp[-\frac{r^2}{4}(9-3\cos 2\theta -\cos 2\phi -3\cos 2\phi\cos 2\theta)]$$
Then using $r=\sqrt{x^2+y^2+z^2}$ and $dxdydz=r^2\sin\phi drd\phi d\theta$ i have rewritten the problem as:
$$V_0=\frac{\rho_0}{4\pi \epsilon_o}\int_0^{\infty}\int_0^{\frac{\pi}{2}}\int_0^{\frac{\pi}{2}}e^{-\frac{9r^2}{4}} e^{-\frac{3r^2}{4}\cos 2\theta} e^{-\frac{r^2}{4}\cos 2\phi} e^{-\frac{3r^2}{4}\cos 2\phi\cos 2\theta}r\sin\phi d\theta d\phi dr$$
Now this just doesn't seem right. I still have multiple variables in most of the  exponents and i don't see how i could separate them. The problem may be that i haven't done any multivariable calculus in over a year but i think i a missing something pretty obvious. Any help would be much appreciated. 
 A: You're completely correct in noticing that spherical coordinates can simplify the equation. The only catch is the method of embedding such a sphere space. Here's what I mean.
Originally, we are given the space:
$$B=\left\{\begin{matrix}x=x,x=[0,\infty] \\ y=y, y=[0,\infty] \\ z=z, z=[0,\infty]\end{matrix}\right.$$
One noticible thing about the parametrization of the sphere you've given:
$$\Omega_1=\left\{\begin{matrix}x=r\sin\phi\cos\theta \\ y=r\sin\phi\sin\theta \\ z=r\cos\phi\end{matrix},\phantom{..}\begin{matrix}r=[0,\infty] \\ \theta=[0..1/2\pi] \\ \phi=[0..1/2\pi]\end{matrix}\right.$$
Is that you'll find that:
$$x^2+y^2=r^2\sin^2\phi\phantom{..}\&\phantom{..}x^2+y^2+z^2=r^2$$
By expanding the integral, we find that:
$$V_0=\frac{\rho_0}{4\pi\epsilon_0}\int_B\frac{\text{e}^{-(x^2+2y^2+2z^2)}}{\sqrt{x^2+y^2+z^2}}dB=\frac{\rho_0}{4\pi\epsilon_0}\int_B\frac{\text{e}^{-(x^2+y^2+z^2)}\text{e}^{-(y^2+z^2)}}{\sqrt{x^2+y^2+z^2}}dB$$
We see that instead of having $x^2+y^2=r%2\sin^2\phi$, we need to have $z^2+y^2=r^2\sin^2\phi$. Thus we reformulate the space to the following parametrization:
$$\Omega_2=\left\{\begin{matrix}x=r\cos\phi \\ y=r\sin\phi\sin\theta \\ z=r\sin\phi\cos\theta\end{matrix},\phantom{..}\begin{matrix}r=[0,\infty] \\ \theta=[0..1/2\pi] \\ \phi=[0..1/2\pi]\end{matrix}\right.$$
Since all the coordinates are equivalent, the angular parametrizations are equivalent. In addition, since this is describing a sphere, the volume element is the same. This produces the following integral:
$$V_0=\frac{\rho_0}{4\pi\epsilon_0}\int_{\Omega_2}\frac{\text{e}^{-r^2}\text{e}^{-r^2\sin^2\phi}}{r}d\Omega_2=\frac{\rho_0}{4\pi\epsilon_0}\int_{\Omega_2}\frac{\text{e}^{-r^2(1+\sin^2\phi)}}{r}d\Omega_2$$
Implementing all limits, the final integral is much simpler than the first as:
$$\frac{\rho_0}{4\pi\epsilon_0}\int^{\frac{1}{2}\pi}_0\int^{\frac{1}{2}\pi}_0\int^\infty_0\frac{\text{e}^{-r^2(1+\sin^2\phi)}}{r}r^2\sin\phi dr d\phi d\theta$$
Implementing all formula, we find that:
$$\frac{\rho_0}{4\pi\epsilon_0}\int^{\frac{1}{2}\pi}_0\int^{\frac{1}{2}\pi}_0\int^\infty_0r\text{e}^{-r^2(1+\sin^2\phi)}\sin\phi dr d\phi d\theta$$
We can also allow $u=r^2$ to finally simplify the integral to:
$$V_0=\frac{\rho_0}{8\pi\epsilon_0}\int^{\frac{1}{2}\pi}_0\int^{\frac{1}{2}\pi}_0\int^\infty_0\text{e}^{-u(1+\sin^2\phi)}\sin\phi du d\phi d\theta$$
Evaluating past the first integral, we find:
$$V_0=\frac{\rho_0}{8\pi\epsilon_0}\int^{\frac{1}{2}\pi}_0\int^{\frac{1}{2}\pi}_0\frac{\sin\phi}{1+\sin^2\phi} d\phi d\theta=\frac{\rho_0}{16\epsilon_0}\int^{\frac{1}{2}\pi}_0\frac{\sin\phi}{1+\sin^2\phi} d\phi=\frac{\rho_0}{16\epsilon_0}\int^{\frac{1}{2}\pi}_0\frac{\sin\phi}{2-\cos^2\phi} d\phi$$
Using the subs $v=\cos\phi$, we finally obtain the single integral:
$$V_0=\frac{\rho_0}{16\epsilon_0}\int_0^1\frac{1}{2-v^2} du$$
I trust you can finish off the integral from here! It was just a matter of finding the right parametrization for the specific charge density
