How can I modify this integral for a numerical integration? I have the following integral$$\int_0^\infty r^k \psi(r) dr \tag{1}$$where $$\psi=e^{-ar-br^2-cr^3-dr^4}\tag{2}$$ 
This is not a well-behaved integrand for numerical integration and diverges in different range of coefficients. So how can I convert it to a well-behaved one? I found Gauss–Laguerre quadrature which uses the following approximation $$\int_0^{\infty}e^{-x}f(x)dx\approx\sum_{i=1}^n w_if(x_i) \tag{3}$$However I'm not sure that this method be usable for my case, because I have 4 terms in exponential which their coefficients maybe positive depending on the initial conditions. So How can address this problem? what is the best way to do numerical integration for (1)?
 I will be so grateful if someone suggests a solution.
 A: Let $I(a,b,d,c,k) = \int_0^\infty \; r^k \mathrm{e}^{-ar -br^2 -cr^3 - dr^4} \,\mathrm{d}r$.


*

*If $d<0$, $I$ diverges.

*If $d = 0$,


*

*If $c < 0$, $I$ diverges.

*If $c = 0$,


*

*If $b < 0$, $I$ diverges.

*If $b = 0$,


*

*If $a < 0$, $I$ diverges.

*If $a = 0$, $I$ diverges.

*If $a > 0$, 


*

*If $k$ is a negative integer, $I$ diverges.

*If $k$ is not a negative integer, $I = a^{1-k} \Gamma(k+1)$.



*If $b > 0$,


*

*If $k$ is a negative integer, $I$ diverges.

*If $k$ is not a negative integer, ...





When $I$ does not diverge, try Gauss-Laguerre quadrature.  This method is designed for semi-infinite intervals of integration and is well-behaved for the $r^k$ in the integrand.  (Also, the Raboniwitz and Weiss paper that describes the method in detail is available.)  It might be helpful to write
$$  \int_0^\infty r^k \mathrm{e}^{-ar-br^2-cr^3-dr^4} \,\mathrm{d}r \\= \int_0^\infty r^k \mathrm{e}^{-r}\mathrm{e}^{(1-a)r-br^2-cr^3-dr^4} \,\mathrm{d}r  \text{.}$$ 
If I were to construct a method from scratch, I would study the polynomial $-dr^4 - cr^3-br^2 - ar$, of degree $n$, to find its maxima.  Then use the estimate $k \ln r < r+1$ to estimate the maxima of $-dr^4 - cr^3-br^2 - ar + k \ln r$ (since this is what we get from $r^k = \mathrm{e}^{k \ln r}$).  Then find an interval of the form $[N,\infty)$ on which the integral is too small to be relevant given our precision/accuracy goals by comparison with $\int_N^\infty r^k \mathrm{e}^{-r^n} \,\mathrm{d}r$.  This integral can be bound/estimated by pretty much any method you like (since the value is $\frac{1}{n} \Gamma\left( \frac{k+1}{n}, N^n \right)$, essentially an incomplete Gamma function evaluation).
If there is more than one maximum, determine whether the minimum between them sits in a valley small enough to contain an ignorable interval, again by comparison with $r^k \mathrm{e}^{-r^n}$.
This leaves potentially two maxima to sample finely enough to reach our precision/accuracy goals and potentially an interval $[0,m]$ on which the function is either small or large, depending on the sign of $k$.  For this interval on the left, if the function is small, sample sparsely enough to reach our precision/accuracy goal.  Otherwise, try Tanh-Sinh quadrature on $[0,m]$.
