The point of this derivation is to expand the exponent, and convert the integral into an infinite sum of integrals, each of which is basically a $\Gamma$-function. Then the sum can be written as a sum of two hypergeometric series. I'm pretty sure this is what Mathematica does.
Call the integral $I$, and write the exponential as
$$ \exp(-(ax-b)^2) = e^{-b^2-a^2x^2}e^{2abx} = e^{-b^2-a^2x^2}\sum_{k\geq0} (2ab)^k\frac{x^k}{k!}. $$
Then, using Mathematica, we have
$$\int_0^\infty e^{-a^2x^2}x^{k+p}dx = \frac{a^{-1-k-p}}{2} \Gamma\left(\frac{1+k+p}{2}\right). $$
After a little manipulation the integral becomes
$$ I = \frac{1}{2}e^{-b^2}a^{-1-p}\sum_{k\geq0} (2b)^k \Gamma\left(\frac{k}{2} + \frac{1+p}{2}\right) \frac{1}{k!}, $$
call the sum $S$.
Now, by definition we have
$$ F(a;b;x) = \sum_{k\geq0}\frac{\Gamma(a+k)/\Gamma(a)}{\Gamma(b+k)/\Gamma(b)}\frac{x^k}{k!}, $$
where $\Gamma(a+k)/\Gamma(a)$ are the rising powers of $a$,
$$ \frac{\Gamma(a+k)}{\Gamma(a)} = a(a+1)(a+2)\cdots(a+k-1). $$
Split $S$ into two sums, one over even $k$, the other over odd $k$; this gets rid of $k/2$ in the argument of $\Gamma$. Then
$$ S = \sum_{k\geq0}\left((2b)^{2k} \frac{\Gamma\left(k+\frac{1+p}{2}\right)}{(2k)!} + (2b)^{2k+1}\frac{\Gamma\left(k+\frac{2+p}{2}\right)}{(2k+1)!}\right). $$
Using
$$ \frac{2^{2k}}{(2k)!} = \frac{1}{k!}\frac{1}{\Gamma(\frac{1}{2}+k)/\Gamma(\frac{1}{2})}, \qquad
\frac{2^{2k}}{(2k+1)!} = \frac{1}{k!}\frac{1}{\Gamma(\frac{3}{2}+k)/\Gamma(\frac{3}{2})}
$$
Find that
$$ S = \Gamma\left(\frac{1+p}{2}\right) F\left(\frac{1+p}{2};\frac12;b^2\right) + 2b\Gamma\left(1+\frac{p}{2}\right)F\left(1+\frac{p}{2};\frac32;b^2\right). $$
Substitute back, use $\Gamma(1+p/2)=(p/2)\Gamma(p/2)$, and you get the required expression for $I$.