To complement Dr. MV's answer I'd like to prove that the differential equation in $f$,
$$
0 = f''(x)-2xf'(x)-2(1+\beta)f(x), \tag{1}
$$
has a solution with an asymptotic expansion of the form
$$
f(x)\sim \sum_{n=0}^\infty a_nx^{-(n+r)} \qquad \text{as } x \to \infty \tag{2}
$$
for some $r \in \mathbb R$. I'll use the same method that is used to calculate the asymptotics of the Airy functions.
To begin, we will search for a solution of $(1)$ of the form
$$
f(x) = \int_\gamma g(z) e^{xz}\,dz, \tag{3}
$$
where $\gamma$ is a (possibly infinite) contour in the complex plane which will be determined.
We calculate
$$
f'(x) = \int_\gamma z g(z) e^{xz}\,dz \qquad \text{and} \qquad f''(x) = \int_\gamma z^2 g(z) e^{xz}\,dz,
$$
and integrating by parts yields
$$
xf'(x) = \int_\gamma [zg(z)]\cdot[xe^{xz}]\,dz = \Bigl[zg(z) e^{xz} \Bigr]_{z=a}^{z=b} - \int_\gamma [g(z) + zg'(z)] e^{xz}\,dz,
$$
where $a$ and $b$ are the (possibly infinite) endpoints of $\gamma$. Substituting these into $(1)$ we get
$$
0 = \int_\gamma \Bigl[ 2zg'(z) + (z^2-2\beta)g(z) \Bigr] e^{xz}\,dz - \Bigl[zg(z) e^{xz} \Bigr]_{z=a}^{z=b}. \tag{4}
$$
The integral vanishes if
$$
0 = 2zg'(z) + (z^2-2\beta)g(z),
$$
and one solution to this equation is
$$
g(z) = z^\beta e^{-z^2/4}, \tag{5}
$$
where the branch cut for $z^\beta$ is taken to be the negative real axis: $\{x : x<0\}$.
With this $g$ the boundary term in $(3)$, namely
$$
\Bigl[zg(z) e^{xz} \Bigr]_{z=a}^{z=b},
$$
vanishes if $\gamma$ begins and ends at $\infty$ in the quadrants defined by $\operatorname{Re} z^2 > 0$, as shown in the image below:

Taking into account the branch cut of $g$ along $\{x : x < 0\}$ we observe that there are exactly three possible contours (up to deformation and changes in orientation):
$$
\begin{align}
&\gamma_1 : -\infty-i\epsilon \longrightarrow (0+) \longrightarrow -\infty + i\epsilon, \\
&\gamma_2 : -\infty+i\epsilon \longrightarrow \infty, \\
&\gamma_3 : -\infty-i\epsilon \longrightarrow \infty,
\end{align}
$$
where $(0+)$ means that the contour passes around the point $z=0$ with positive orientation. These three contours are shown in the figure below, and the branch cut of $g$ is shown as a black jagged line:

Since the differential equation $(1)$ is of second order we expect that the solutions $f$ given by one of these contours must be dependent on the other two, and indeed we have
$$
\int_{\gamma_3} = \int_{\gamma_1} + \int_{\gamma_2}.
$$
So we can throw the contour $\gamma_3$ out, and we are finally left with two linearly independent solutions of $(1)$ which we call $f = F_1$ and $f = F_2$, given by
$$
\begin{align}
F_1(x) &= \int_{\gamma_1} z^\beta \exp(xz - z^2/4)\,dz, \\
F_2(x) &= \int_{\gamma_2} z^\beta \exp(xz - z^2/4)\,dz.
\end{align} \tag{6}
$$
I'll omit the analysis of $F_2(x)$ as $x \to \infty$, but suffice it to say that a relatively straightforward application of the saddle point method yields the asymptotic expansion
$$
F_2(x) \sim (2x)^\beta \exp(x^2) \sum_{n=0}^{\infty} \binom{\beta}{2n} \Gamma\!\left(n + \frac{1}{2}\right) x^{-2n} \qquad \text{as } x \to \infty.
$$
To leading order,
$$
F_2(x) \sim \sqrt{\pi} (2x)^\beta e^{x^2} \qquad \text{as } x \to \infty.
$$
Since $F_2(\infty) \neq 0$ this is not the solution we're looking for. Let's take a closer look at $F_1$.
Assuming $x > 0$ we calculate
$$
\begin{align}
F_1(x) &= \int_{\gamma_1} z^\beta \exp\!\left(xz - z^2/4\right) dz \\
&= e^{-i\pi \beta} \int_{-\infty}^0 (-z)^\beta \exp\!\left(xz - z^2/4\right) dz - e^{i\pi \beta} \int_{-\infty}^0 (-z)^\beta \exp\!\left(xz - z^2/4\right) dz \\
&= -2i\sin(\beta \pi) \int_0^\infty t^\beta \exp\!\left(-xt - t^2/4\right) dt \\
&= -2i\sin(\beta \pi) x^{\beta+1} \int_0^\infty \exp\!\left[-x^2 \left(u+u^2/4\right)\right] du \\
&= -i\sin(\beta \pi) (2x)^{\beta+1} \int_0^\infty \frac{\left(\sqrt{1+v} - 1\right)^\beta}{\sqrt{1+v}} \exp\!\left(-x^2v\right) dv,
\end{align}
$$
where we have made the successive substitutions $z = -t$, $t=xu$, and $u+u^2/4 = v$. The non-exponential factor in the integrand can be represented by a convergent power series,
$$
\frac{\left(\sqrt{1+v} - 1\right)^\beta}{\sqrt{1+v}} = v^\beta \sum_{n=0}^{\infty} a_n v^n, \qquad |v| < 1,
$$
so an application of Watson's lemma tells us that $F_1(x)$ has an asymptotic expansion of the form
$$
F_1(x) \sim x^{-\beta+1} \sum_{n=1}^{\infty} b_n x^{-2n} \qquad \text{as } x \to \infty
$$
for some coefficients $b_n$. To leading order,
$$
F_1(x) \sim B x^{-\beta-1} \qquad \text{as } x \to \infty
$$
for some constant $B$.
We have therefore shown that the differential equation $(1)$ has a solution with an asymptotic expansion of the form $(2)$, and in principle we could calculate the coefficients from this integral representation of the solution.