Im looking at a condensate in a spherically symeric trap so the potential is $$V(r)=\frac{m}{2}\omega_r^2r^2$$ When the interactions are weak we assume this wavefunction:
$$\psi(\boldsymbol{r})=\left(\frac{N}{\pi^{\frac{3}{2}}\sigma^3l_r^3}\right)^{\frac{1}{2}}\exp\left(-\frac{r^2}{2\sigma^2l_r^2}\right)$$
with $r^2=x^2+y^2+z^2$ and $l_r=\sqrt{\frac{\hbar}{m\omega_r}}$ the energy is given by $$E=\int\frac{\hbar}{2m}|\nabla\psi|^2+V|\psi|^2+\frac{g}{2}|\psi|^4 \mathrm{d}^3\boldsymbol{r}$$ with $g=4\pi\hbar^2a_s/m$ I know the solution should be
$$E(\sigma)=\hbar\omega_rN\left[\frac{3}{4\sigma^2}+\frac{3\sigma^2}{4}+\frac{1}{\sqrt{2\pi}}\frac{Na_s}{l_r}\frac{1}{\sigma^3}\right]$$
The problem is the exponential term that pops up in each term of the energy integral casue you end up with a bunch of error functions and it doesnt cancel down nicely is there something im missing?
Apologies, I have been busy the last couple days. For brevity's sake I will rewrite $$\psi(r)=\alpha\cdot \exp(-\beta r^2)$$ We can simplify the problem a little bit by only looking at the one dimensional integral (assuming the bounds are $[0,\infty)$ as per usual in quantum mechanics) $$\int_0^\infty ...~\mathrm{d}r$$ Since $\psi$ is radially symmetric, we can take this as a constant and then just integrate that over $\theta\in[0,2\pi)$, $\phi\in[0,\pi)$. Recall that if $\psi$ has no $\theta,\phi$ dependence, we can write simply $$\nabla \psi=\frac{\mathrm{d} \psi}{\mathrm{d} r}=-2\alpha\beta r\exp(-\beta r^2)$$ Which leads us to the one dimensional version of the integral in the question $$E_1=\int_0^\infty \frac{\hbar}{m}2\alpha^2\beta^2r^2\exp(-2\beta r^2)+\frac{m\omega^2}{2}r^2\alpha^2\exp(-2\beta r^2)+\frac{g}{2}\alpha^4\exp(-4\beta r^2)\mathrm{d}r$$ If you're prepared to accept the integral identity $$\int_0^\infty x^n\exp(-kx^2)\mathrm{d}x=\frac{1}{k^{(n+1)/2}}\frac{1}{2}\Gamma\left(\frac{n+1}{2}\right)$$ Which is discussed in this MSE question, then we can write this as $$E_1=\frac{1}{(2\beta)^{3/2}}\frac{1}{2}\Gamma(3/2)\left(\frac{\hbar}{m}2\alpha^2\beta^2+\frac{m\omega^2\alpha^2}{2}\right)+\frac{g\alpha^4}{2}\frac{1}{(4\beta)^{1/2}}\frac{1}{2}\Gamma(1/2)$$ The half integer values of the Gamma are well known. See Wikipedia. So we can make some simplifications $$E_1=\frac{\sqrt{\pi}\alpha^2}{2^{5/2}}\left(\frac{\hbar\sqrt{\beta}}{m}+\frac{m\omega^2}{4\beta^{3/2}}+\frac{g\alpha^2}{\sqrt{2}\sqrt{\beta}}\right).$$