How to derive the electron density of Hooke's atom?

[Hooke's atom](https://en.wikipedia.org/wiki/Hooke%27s_atom) has the following wave function:


$$\Psi( \boldsymbol{r\_1}, \boldsymbol{r\_2}) \propto \left( 1 + \frac{|\boldsymbol{r\_1} - \boldsymbol{r\_2}|}{2} \right) \exp\left(-\frac{(r\_1^2+r\_2^2)}{4}\right).$$


Question: How to calculate electronic density $\rho (\boldsymbol{r\_1})$ using this wave function?


The answer, which should contain Error function, is on [Wikipedia](https://en.wikipedia.org/wiki/Hooke%27s_atom) but not method.


I tried $\rho (\boldsymbol{r\_1}) = \int |\Psi (( \boldsymbol{r\_1}, \boldsymbol{r\_2}))|^2 d\boldsymbol{r\_2}$ but got wrong answer. In spherical coordinates:


$\rho (\boldsymbol{r\_1}) = \int |\Psi (( \boldsymbol{r\_1}, \boldsymbol{r\_2}))|^2 d\boldsymbol{r\_2} = \int\_0^{2 \pi} \int\_0^{\pi} \int\_0^{\infty} |\Psi (( \boldsymbol{r\_1}, \boldsymbol{r\_2}))|^2 r\_2^2 d{r\_2} \sin(\theta\_2) d\theta\_2 d\phi\_2$


$\rho (\boldsymbol{r\_1}) \propto \int\_0^{2 \pi} \int\_0^{\pi} \int\_0^{\infty} \left( 1 + \tfrac{1}{2} \sqrt{r\_1^2 + r\_2^2 -2 \boldsymbol{r\_1} \cdot \boldsymbol{r\_2} } \right)^2 \exp\left(-\frac{(r\_1^2+r\_2^2)}{2}\right) r\_2^2 d{r\_2} \sin(\theta\_2) d\theta\_2 d\phi\_2$


where $\boldsymbol{r\_1} \cdot \boldsymbol{r\_2} = r\_1 r\_2 \left[ \sin(\theta\_1) \sin(\theta\_2) \cos(\phi\_1) \cos(\phi\_2) + \sin(\theta\_1) \sin(\theta\_2) \sin(\phi\_1) \sin(\phi\_2) + \cos(\theta\_1) \cos(\theta\_2) \right]$


Aso according to the article linked on wikipedia answer should contain Error function. However Mathematica does not yield this result upon carrying out above integral. Where is my mistake?


(Note: crosspost from [physics.se](https://physics.stackexchange.com/questions/284350/electronic-density-of-hooke-atom)).


