Integral evaluation through series expansion

$$I(\gamma,b,\beta)=\frac{\beta}{\sqrt{\beta^2+b^2}}\int_{0}^{\infty} \frac{1}{x^2+\gamma^2\,\frac{\beta^2}{b^2+\beta^2}} \, e^{- \sqrt{\beta^2+b^2} \, \sqrt{x^2+\gamma^2}} \, \mathrm{d}x$$

I was wondering if it makes any sense to use the power series of

$$e^{- \sqrt{\beta^2+b^2} \, \sqrt{x^2+\gamma^2}} = \sum_{n=0}^\infty \frac{(- \sqrt{\beta^2+b^2} \, \sqrt{x^2+\gamma^2})^n}{n!}$$

and the indefinite integral (from Mathematica)

$$I_1(x,\gamma,b,\beta,n)=\frac{\beta}{\sqrt{\beta^2+b^2}}\int \frac{(- \sqrt{\beta^2+b^2} \, \sqrt{x^2+\gamma^2})^n}{\left(x^2+\gamma^2\,\frac{\beta^2}{b^2+\beta^2}\right)n!}\, \mathrm{d}x=\frac{(-1)^n}{\beta\,\gamma^2\,n!}\,x\,(\beta^2+b^2)^{\frac{1+n}{2}}\,\left(1+\frac{x^2}{\gamma^2}\right)^{-n/2}\,(x^2+\gamma^2)^{\frac{n}{2}}\,\mathrm{F}_1\left(\tfrac{1}{2};-\tfrac{n}{2},1,\tfrac{3}{2};-\tfrac{x^2}{\gamma^2},-\tfrac{x^2\,(b^2+\beta^2)}{\beta^2\,\gamma^2}\right)$$

with $\mathrm{F}_1$ the Appell Hypergemetric Function.

Since $$\lim_{x\to0} I_1(x,\gamma,b,\beta,n)=0$$ I would guess that $$I(\gamma,b,\beta)= \lim_{x\to\infty}\sum_{n=0}^\infty \frac{(-1)^n}{\beta\,\gamma^2\,n!}\,x\,(\beta^2+b^2)^{\frac{1+n}{2}}\,\left(1+\frac{x^2}{\gamma^2}\right)^{-n/2}\,(x^2+\gamma^2)^{\frac{n}{2}}\,\mathrm{F}_1\left(\tfrac{1}{2};-\tfrac{n}{2},1,\tfrac{3}{2};-\tfrac{x^2}{\gamma^2},-\tfrac{x^2\,(b^2+\beta^2)}{\beta^2\,\gamma^2}\right)$$

If this is even correct, is there a way to get a closed form solution from this approach? I know there is a solution for $b=0$. Thank you.

Edit: I think it can be simplified to $$I(\gamma,b,\beta)= \lim_{x\to\infty}\sum_{n=0}^\infty \frac{(-1)^n}{\beta\,\gamma^{2-n}\,n!}\,x\,(\beta^2+b^2)^{\frac{1+n}{2}}\,\,\mathrm{F}_1\left(\tfrac{1}{2};-\tfrac{n}{2},1,\tfrac{3}{2};-\tfrac{x^2}{\gamma^2},-\tfrac{x^2\,(b^2+\beta^2)}{\beta^2\,\gamma^2}\right)$$

First, let's get the parameters in a more convenient form:

$$I=\frac{\beta}{p} \int_0^\infty \frac{e^{-p \sqrt{t^2+1}}}{t^2+q^2}dt$$

Where:

$$p=\gamma \sqrt{\beta^2+b^2}, \qquad q^2=\frac{\beta^2}{b^2+\beta^2} \leq 1$$

Now consider the new integral, and let's make another substitution:

$$t=\sinh u$$

$$\int_0^\infty \frac{e^{-p \sqrt{t^2+1}}}{t^2+q^2}dt=\int_0^\infty \frac{e^{-p \cosh u}\cosh u}{\cosh^2 u-s^2}du=f(p,s)$$

Where:

$$s^2=1-q^2 \leq 1$$

We can expand the denominator as a geometric series for the full range of the variable, getting the Taylor series in $s$:

$$f(p,s)= \sum_{k=0}^\infty s^{2k} \int_0^\infty \frac{e^{-p \cosh u}}{\cosh^{2k+1} u}du$$

The integrals in general don't have a closed form (though they all converge, so can be computed numerically). However, we can obtain a differential equation in $p$.

Let's now denote $f(p,s)=g(p)$, leaving $s$ as a general parameter. Then:

$$g''(p)=\sum_{k=0}^\infty s^{2k} \int_0^\infty \frac{e^{-p \cosh u}}{\cosh^{2k-1} u}du= \\ =\int_0^\infty e^{-p \cosh u} \cosh u~du+s^2 \sum_{k=0}^\infty s^{2k} \int_0^\infty \frac{e^{-p \cosh u}}{\cosh^{2k+1} u}du$$

The first integral has a closed form in terms of a modified Bessel function. Finally we obtain a second order ODE:

$$g''(p)-s^2 g(p)-K_1 (p)=0 \tag{1}$$

The initial conditions are a little more complicated. From the original integral we have:

$$g(0)=\frac{\pi}{2 \sqrt{1-s^2}}$$

We can use $g( \infty)=0$ as another condition, however it's not suitable for numerical solution of the ODE. In that case, we can just numerically compute some value, like $g(1)$ or $g'(1)$. The derivative at $0$ diverges.

The equation has an analytic solution, but not in closed form. The homogeneous part obviously gives an exponential solution, but the Bessel function makes the final solution much more complicated.

Here's a numerical solution of the ODE (blue) compared to numerical integration (orange). Both are done in Mathematica for $s=e^{-1}$:

For $s=1$, Wolfram Alpha and Mathematica give a closed form general solution to the ODE in herms of Meijer G-function:

http://www.wolframalpha.com/input/?i=g%27%27%5Bx%5D-g%5Bx%5D-BesselK%5B1,x%5D%3D0

$$g(p)=-\frac{1}{4 \sqrt{\pi }} \left( \pi e^{p} G_{2,3}^{3,0}\left(2 p\left| \begin{array}{c} 1,\frac{3}{2} \\ 0,0,2 \\ \end{array} \right.\right)+e^{-p}G_{2,3}^{3,1}\left(2 p\left| \begin{array}{c} \frac{3}{2},1 \\ 0,0,2 \\ \end{array} \right.\right)\right)+c_1 e^p+c_2 e^{-p}$$

$$s=1$$

Note that for $s=1$ the original integral diverges.

As such, I doubt there's anything more simple, especially for general $s$.