I think the main problem for me is to calculate the integral of $$\int_{0}^{\infty}\frac{e^{-\frac{r^2}{2t}}}{\sqrt{x^2+r^2}}r^{n-1}dr,n\geq2$$
For n=2, change of variable $y=\sqrt{x^2+r^2}$ would work well. But I'm stuck when n>2.
The background of this problem is to calculate the expectation of n+1-dimension Inverse Bessel Process starts at $(1,0,...0)$, where the orginal integral is $$(2\pi t)^{-(n+1)/2}\int_{\mathbb{R}^{n+1}}\frac{e^{-||\overrightarrow{x}-(1,0,...0)||^2/(2t)}}{||\overrightarrow{x}||}dx$$
I try to first integrate the last $n$ coordinates and change them into spherical coordinates, taking out the constant, if I'm correct, you will get the integral at the beginning.
You are free to use functions like erf() or cdf of Normal distribution or Gamma function to express your result, if no such expression exists, a recursive formula is also desired.