specific improper integral involving erf 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$).
 A: Here's some ideas. Your integrand is non-oscillatory with good decay so you should only need <100 quadrature pts for machine accuracy if you use a high-order method. I suggest plain trapezoid rule with Alpert end correction at the bottom end and truncation in the Gaussian tail at top end. Eg translate code from http://mpspack.googlecode.com/svn/trunk/@quadr/QuadNodesInterval.m
http://mpspack.googlecode.com/svn/trunk/@quadr/QuadSmoothExtraPtNodes.m
If erf is slow build a look-up table using piecewise polynomial interpolation on Chebyshev grids.
Also: sub y=px and your function of $x$ is an antiderivative of $(\mbox{erf} y)^\alpha e^{-y^2} / y$; if this is to be summed over $x$ as you say, compute the antiderivative together as you go along the $x$ integral. 
Also: try to find recurrence relation for $\alpha$ to $\alpha+1$ if you care about many alpha cases (I didn't spot anything yet).
A: As others have pointed out, there's no hope for a closed-form solution. To evaluate it numerically in an efficient way, I would first rewrite this as
$$
\int_x^\infty (\mathrm{erf}(y))^\alpha \exp(-y^2) {dy \over y}\;.
$$
Now if $y$ is larger than about $10$ or so, then the integrand is equal to $0$ for all practical purposes, so that the integral is pretty well approximated by
$$
\int_x^{10} (\mathrm{erf}(y))^\alpha \exp(-y^2) {dy \over y}\;,
$$
which can easily be evaluated to machine precision by any quadrature formula (provided you can evaluate the error function). If you want your code to be fast, just precompute the integrand at a few thousand grid points...
If you are interested in small relative error rather than just small absolute error, then you can also get asymptotics for large $x$. To leading order, you'll get a pretty good approximation by 
$$
\int_x^\infty (\mathrm{erf}(y))^\alpha \exp(-y^2) {dy \over y}
\approx 
{1\over x} \int_x^\infty \exp(-y^2) dy = {\sqrt{\pi} \, \mathrm{erfc}(x) \over 2x} \approx {e^{-x^2} \over 2x^2}\;.
$$
