I have encountered an integral, and kindly ask for help with a solution. It is beyond my own capabilities, and neither Maple nor Mathematica were of any help: $$ \int_{1}^{\infty} \left[\mathrm{erf}\left(px\right)\right]^{\alpha} \exp\left(-p^2x^2\right)\,\frac{\mathrm{d}p}{p}, $$ where $\alpha$ is an integer $> 0$, and $x>0$.

If the lower integration boundary were $0$ instead of $1$, a Laplace transform might be of help (with the transformation $p\rightarrow k=p^2x^2$). However, integrating that part separately leaves me with another integral I can't tackle (with integration boundaries $0$ and $x^2$, or $0$ and $1$ before the transformation).

I would also be grateful for hints for efficiently solving this integral numerically if it cannot be simplified further. Again, for $0$ as the lower integration boundary, (generalized) Gauss-Laguerre quadrature would help for one part, but the other is still a problem (presumably Gauss-Legendre quadrature is not conducive because of the singularity at $0$).