Your particular result is true because it happens to be a Laplace Transform of a function that does not grow too quickly near the postiive real line.
Indeed, your book of Gamelin (page 367) expresses the trigamma $\psi'$ in this form,
$$ \psi'(z)=\frac{d}{dz} \frac{\Gamma'(z)}{\Gamma(z)} = \int_0^\infty e^{-zt}\frac{t}{1-e^{-t}} dt, \quad \Re z> 0.$$
This works because $$f(t) = \frac{t}{1-e^{-t}}$$ is exponential type of order $0+$ on a small neighbourhood of the positive real line, which fits into a more general result:
Theorem Let $\delta>0$, $\rho,\rho'>0$, with $\rho'>\rho$.
Suppose $U_\delta$ is a $\delta$-neighbourhood of $\mathbb R_+=(0,\infty)$ in $\mathbb C$, i.e. $$ U_\delta = \{t\in\mathbb C : d(t,\mathbb R_+) < \delta \}.$$ If $f$ is analytic of exponential order $\rho$ in $U_\delta$, then
$$ \mathcal Lf(z) = \int_0^\infty e^{-zt}f(t) dt \simeq \sum_{k=0}^\infty \frac{f^{(k)}(0)}{z^{k+1}} $$
as $z\to \infty$ in the set $\{\Re z > \rho'\}$.
Implicit above is the fact that $\int_0^\infty e^{-zt}f(t) dt$ converges for analytic functions of exponential order $\rho$ in $U_\delta$. We say $f$ is of exponential order $\rho>0$ in $D$ if there exists $M$ such that for every $t\in D$,
$$ |f(t)| \le M e^{\rho |t|}.$$
An example: if $f(t)=1$ then we can take any $\rho>0$ and $\mathcal Lf = \frac{1}{z}\simeq \frac1z + 0 + \dots $ is indeed valid on every half plane $ \{\Re z > \rho' \}$.
The proof is relatively straightforward; here's the sketch. The appearance of derivatives strongly suggests integration by parts. Under the given assumptions, you can check that $f^{(k)}$ is also of exponential type $\rho$ for every $k$, which justifies the integration by parts. We record the constant that verifies this as $M_k$, i.e.
$$ |f^{(k)}(z)| \le M_k e^{\rho |z|}.$$
Performing the integration by parts inductively leads to
$$\mathcal{L} f(z)=\sum_{k=0}^{n-1} \frac{f^{(k)}(0)}{z^{k+1}}+\frac{1}{z^n} \int_{0}^{\infty} e^{-z t} f^{(n)}(t) d t$$
So we need to bound this last integral. One last integration by parts gives
\begin{align}
\left| \int_{0}^{\infty} e^{-z t} f^{(n)}(t) d t\right|
& \le \frac{|f^{(n)}(0)|}{|z|} + \frac1{|z|}\int_0^\infty \underbrace{|e^{-tz} f^{(n+1)} (t)|}_{\Large\le M_{n+1} e^{-t \Re z+\rho t}} dt
\\
&\le \frac{M_n}{|z|} + \frac{M_{n+1}}{|z|(\Re z-\rho)}
\\
&\le \frac{M_n}{|z|} + \frac{M_{n+1}}{|z|(\rho'-\rho)}
\\&=: \frac{C_{n+1}}{|z|}
\end{align}
where $\rho<\rho'<\Re z$ was crucially used. This gives the required estimate to satisfy the definition of an asymptotic series (with asymptotic sequence $(z^{-k-1})_{k\ge 0}$).
Here is the bibliography from the description of my old course on Asymptotic Analysis, may be useful...
Books:
We will not follow any particular book, but most of the material can be found in:
- C.F. Carrier, M. Krook and C.E. Pearson, Functions of a Complex Variable: theory and technique, Hodbooks.
- N.G. De Bruijn, Asymptotic Methods in Analysis, North-Holland Publishing co. (3d ed.) (1970).
- P.P.G. Dyke, An Introduction to Laplace Transforms and Fourier Series, Springer Undergraduate Mathematics Series (2000).
- G. Hardy, Divergent Series, Clarendon Press, 1963/American Mathematical Society, 2000.
- R.B. Dingle, Asymptotic Expansions: Their Derivation and Interpretation, Academic Press (1973).
So as mentioned in the comments this identity provides a "closed form" for your hypergeometric function in terms of Laguerre polynomials. Here is another option that may be considered a more elementary form:
Using DLMF 13.3.17 for $n=2m+1$ and $a=1$ you may write
$$
\begin{aligned}
{_1F}_1\left({2m+2\atop 1/2};z\right)%
&=\frac{1}{(2m+1)!}\partial_z^{2m+1} z^{2m+1} {_1F}_1\left({1\atop 1/2};z\right)\\
&=\frac{1}{(2m+1)!}\partial_z^{2m+1} z^{2m+1} (1+\sqrt{\pi z} e^z\operatorname{erf}(\sqrt z))\\
&=\frac{1}{(2m+1)!}\left((2m+1)!+\partial_z^{2m+1} z^{2m+1}\sqrt{\pi z} e^z\operatorname{erf}(\sqrt z)\right)\\
&=1+\frac{1}{(2m+1)!}\partial_z^{2m+1} z^{2m+1}\sqrt{\pi z} e^z\operatorname{erf}(\sqrt z)\\
\end{aligned}
$$
You may then apply the General Leibniz rule to evaluate the derivative. Once the derivatives are evaluated you may then make the substitution $z\mapsto-x^2$ to obtain the final result.
There are many other different forms you can put your hypergeometric function into using similar methods and variants of this method. What its going to come down to is what form is useful for your application.
Best Answer
From $\operatorname{erf}(z)=1-\operatorname{erfc}(z)$ and $\S7.12(\text{i})$ in the DLMF, it follows that $$ \operatorname{erf}(z) \sim 1 - \frac{{\mathrm{e}^{ - z^2 } }}{{\sqrt \pi }}\sum\limits_{m = 0}^\infty {( - 1)^m \frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} $$ if $z\to \infty$ and $\operatorname{Re}(z)>0$, and $$ \operatorname{erf}(z) \sim - 1 - \frac{{\mathrm{e}^{ - z^2 } }}{{\sqrt \pi }}\sum\limits_{m = 0}^\infty {( - 1)^m \frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} $$ if $z\to \infty$ and $\operatorname{Re}(z)<0$. The positive and negative imaginary axes are Stokes lines. On the Stokes lines, the correct expansion is the average of the expansions on the two sides, i.e., $$ \operatorname{erf}(z) \sim - \frac{{\mathrm{e}^{ - z^2 } }}{{\sqrt \pi }}\sum\limits_{m = 0}^\infty {( - 1)^m \frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} $$ if $z\to \infty$ and $\operatorname{Re}(z)=0$. Since $\operatorname{erfi}(z)=-\mathrm{i}\operatorname{erf}(\mathrm{i}z)$, we can infer that $$ \operatorname{erfi}(z) \sim \operatorname{sgn}(\operatorname{Im}(z))\mathrm{i} + \frac{{\mathrm{e}^{z^2 } }}{{\sqrt \pi }}\sum\limits_{m = 0}^\infty {\frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} $$ as $z\to \infty$, with the usual convention that $\operatorname{sgn}(0)=0$.