I wonder if there is a quick numerical integration method on the integral below:
$$ \beta = \int_c^\infty (y+b)^{-p}\exp\left(-y^2\right)\ \mathrm{d}y, $$
where $b$ and $c$ are real positive and $p$ is a positive integer.
I have tried decomposing the $(y+b)^{-p}$ term into its Taylor series around $y=0$ and finish the integral. I obtained
$$ \beta = \frac{1}{2} \sum_{n=0}^{\infty} {{p+n-1}\choose n} b^{-(p+n)} (-1)^n \Gamma\left[\frac{1+n}{2}, c^2\right] $$
where $\Gamma[a,z]$ is the upper incomplete gamma function. The problem with the expression above is that the series diverge for small $b$.
Does anyone have an idea on getting a computable form of the integral? At the moment, numerical quadrature method is not an option for me (e.g. trapezoidal, Simpson's, Gaussian, etc).