Integrating over polar coordinates and 6 dimensions I have a wave function as follows:
$$\psi=A \exp (\frac{r}{2a}-br^2-4bR^2) \tag{1}$$where $A$ is normalization constant. Normalization means 
$$A^2\int \psi^2d\tau=1 \tag{2}$$ so $A$ is $$A=\frac{1}{\sqrt{\int \psi^2d\tau}} \tag{3}$$
where $d\tau$ means whole space and in polar coordinates is $$\int_0^{\infty}\int_0^{\pi}\int_0^{2\pi}d\phi\;\sin{\theta}\;d\theta\; r^2dr\tag{4}$$
when I try to obtain $A$ in (1) using (2-4) I get this
$$\frac{\sqrt{\frac{b \exp(-\frac{1}{8a^2b})}{1+Erf(\frac{1}{2a\sqrt{2b}})}}}{\pi^{3/2}} \tag{5}$$ where I have obtained $16\pi^2$ for integrating over  $d\theta_1,d\theta_2,d\phi_1,d\phi_2 $ because we have no arguments in (1) for them and then I have integrated over $r$ and $R$ using Mathematica. But the final result for $A$ in the paper which I try to reproduce its result is $$8(\frac{b}{\pi})^{3/2}\{(\frac{2}{x})\pi^{-1/2}+(1+\frac{2}{x^2})\exp(\frac{1}{x^2})[1+erf(x)]\}^{-1/2} \tag{6}$$ where $$x=2a {(2b)}^{1/2}$$ It's obvious that the general form of my $A$ is correct, but the coefficients don't match. I really don't know how to address this problem! Any idea?
Addendum (For comment 1)
For integrating over $r$ we get
$$\frac{\exp(\frac{1}{8a^2b}-8bR^2)\sqrt{\frac{\pi}{2}}(1+Erf[\frac{1}{2a\sqrt{2b}}])}{2\sqrt{b}} \tag{7}$$ and after integrating (7) over $R$ we obtain
$$\frac{\exp(\frac{1}{8a^2b})\pi^3(1+Erf[\frac{1}{2a\sqrt{2b}}])}{b} \tag{8}$$ and after inversion and get sqrt we reach to final result (5).
Addendum 2
When I evaluate the following expression in Mathematica
$$\int_0^{\infty} r^2 \exp (\frac{r}{a}-2br^2) \tag{9}$$ I get
$$\frac{4a\sqrt{b}+(1+4a^2 b)\exp (\frac{1}{8a^2 b})\sqrt{2\pi}(1+Erf[\frac{1}{2a\sqrt{2b}}])}{64a^2b^{\frac 5 2}} \tag{10}$$
 A: First they worked out
$$\begin{align}I_{\vec R}&=\int_{-\infty}^{\infty}dX\int_{-\infty}^{\infty}dY\int_{-\infty}^{\infty}dZe^{-8b(X^2+Y^2+Z^2)}\\
&=\left(\int_{-\infty}^{\infty}e^{-8bX^2}dX\right)^3=\left(\frac1{\sqrt{8b}}\cdot2\cdot\frac12\operatorname{\Gamma}\left(\frac12\right)\right)^3=\left(\frac{\pi}{8b}\right)^{3/2}\end{align}$$
Then they said
$$\begin{align}I_{\vec r}&=\int_0^{\infty}r^2\,dr\int_0^{\pi}\sin\theta\,d\theta\int_0^{2\pi}d\phi e^{\frac ra-2br^2}=4\pi\left(\frac1{2b}\right)^{3/2}\int_0^{\infty}r^2e^{\frac r{a\sqrt{2b}}-r^2}dr\tag{1}\\
&=4\pi\left(\frac1{2b}\right)^{3/2}e^{\frac1{8a^2b}}\int_0^{\infty}r^2e^{-\left(r-\frac1{2a\sqrt{2b}}\right)^2}dr\tag{2}\\
&=4\pi\left(\frac1{2b}\right)^{3/2}e^{\frac1{8a^2b}}\int_{-\frac1{2a\sqrt{2b}}}^{\infty}\left(r+\frac1{2a\sqrt{2b}}\right)^2e^{-r^2}dr\tag{3}\\
&=4\pi\left(\frac1{2b}\right)^{3/2}e^{\frac1{8a^2b}}\left\{\left[-\frac12\left(r+\frac1{a\sqrt{2b}}\right)e^{-r^2}\right]_{-\frac1{2a\sqrt{2b}}}^{\infty}\right.\tag{4}\\
&\quad\left.+\left(\frac12+\frac1{8a^2b}\right)\int_{-\frac1{2a\sqrt{2b}}}^{\infty}e^{-r^2}dr\right\}\\
&=4\pi\left(\frac1{2b}\right)^{3/2}e^{\frac1{8a^2b}}\left\{\frac1{4a\sqrt{2b}}e^{-\frac1{8a^2b}}\right.\tag{5}\\
&\quad\left.+\left(\frac12+\frac1{8a^2b}\right)\left(\frac{\sqrt{\pi}}2+\frac{\sqrt{\pi}}2\operatorname{erf\left(\frac1{2a\sqrt{2b}}\right)}\right)\right\}\end{align}$$
In $(1)$ I said $\int_0^{\pi}\sin\theta\,d\theta=\left.-\cos\theta\right|_0^{\pi}=2$ and $\int_0^{2\pi}d\phi=\left.\phi\right|_0^{2\pi}=2\pi$. Then I let $r=\frac u{\sqrt{2b}}$.
In $(2)$ I completed the square: $-r^2+\frac r{a\sqrt{2b}}=-\left(r-\frac1{2a\sqrt{2b}}\right)^2+\frac1{8a^2b}$.
In $(3)$ I let $r=u+\frac1{2a\sqrt{2b}}$.
In $(4)$ I integrated $\int_{-1/x}^{\infty}re^{-r^2}dr=\left.-\frac12e^{-r^2}\right|_{-1/x}^{\infty}$ and integrated by parts $\int_{-1/x}^{\infty}r^2e^{-r^2}dr=\left.-\frac12re^{-r^2}\right|_{-1/x}^{\infty}+\frac12\int_{-1/x}^{\infty}e^{-r^2}dr$.
In $(5)$ I said $\int_{-1/x}^{\infty}e^{-r^2}dr=\frac{\sqrt{\pi}}2\left(\frac2{\sqrt{\pi}}\int_0^{\infty}e^{-r^2}dr+\frac2{\sqrt{\pi}}\int_{-1/x}^0e^{-2r^2}dr\right)=\frac{\sqrt{\pi}}2\left(1+\operatorname{erf}\left(\frac1x\right)\right)$.  
With $x=2a\sqrt{2b},$ I get
$$\int\lvert\psi\rvert^2d\tau=I_{\vec R}I_{\vec r}=\frac{\pi^3}{64b^3}\left\{\frac2{x\sqrt{\pi}}+e^{1/x^2}\left(1+\frac2{x^2}\right)\left(1+\operatorname{erf}\left(\frac1x\right)\right)\right\}$$
Is there a typo in your expression where you have $\operatorname{erf}(x)$? OK, so it's a typo in the paper. These kinds of things are hard to get right in a vast sea of typeset equations.
