Integral involving a confluent hypergeometric function I have the following integral involving a confluent hypergeometric function:
$$\int_{0}^{\infty}x^3e^{-ax^2}{}_1F_1(1+n,1,bx^2)dx$$
where $a>b>0$ are real constants, and $n\geq 0$ is an integer.
Wolfram Mathematica returns the following solution: $\frac{a^{n-1}(a+bn)}{(a-b)^{n+2}}$.  However, I can't figure out how it arrived at it (I always try to check the solutions "on paper" that Mathematica gives me -- or at least using Gradshteyn and Ryzhik).  Can anyone help?

Let's start with the hypergeometric function. We have:
\begin{eqnarray}
F_{2,1}[1+n,1;b x^2] &=& 
\sum\limits_{m=0}^\infty  \frac{(1+n)^{(m)}}{m!} \cdot \frac{(b x^2)^m}{m!} \\
&=& \sum\limits_{m=0}^\infty \frac{(m+1)^{(n)}}{n!} \cdot \frac{(b x^2)^m}{m!} \\
&=&\left. \frac{1}{n!} \frac{d^n}{d t^n}  \left( t^n \cdot e^{b x^2 \cdot t} \right)\right|_{t=1} \\
&=& \frac{1}{n!} \sum\limits_{p=0}^n \binom{n}{p} n_{(p)} (b x^2) ^{n-p} \cdot e^{b x^2}
\end{eqnarray}
Therefore the integral in question reads:
\begin{eqnarray}
&&\int\limits_0^\infty x^3 e^{-a x^2} F_{2,1}[1+n,1;b x^2] dx =\\ 
&& \frac{1}{2} \sum\limits_{p=0}^n \frac{n^{(p)}}{p! (n-p)!} b^{n-p} \frac{(n-p+1)!}{(a-b)^{n-p+2}} \\
&&\frac{1}{2} \frac{1}{(a-b)^{n+2}} \sum\limits_{p=0}^n \binom{n}{p} \frac{(n-p+1)^{(p)}}{(n-p+2)^{(p-1)}} \cdot b^{n-p} (a-b)^p \\
&& \frac{1}{2} \frac{1}{(a-b)^{n+2}} \sum\limits_{p=0}^n \binom{n}{p} (n-p+1) \cdot b^{n-p} (a-b)^p \\
&& \frac{1}{2} \frac{1}{(a-b)^{n+2}} a^{n-1} (a+b n)
\end{eqnarray}
