How to find the integral $\int_{0}^{\infty}\exp(- (ax+b/x))\,dx$?

Sub $u=a x+b/x$. Then $a x^2-u x+b = 0$, and therefore

$$x = \frac{u}{2 a} \pm \frac{\sqrt{u^2-4 a b}}{2 a}$$

$$dx = \frac1{2 a}\left (1 \pm \frac{u}{\sqrt{u^2-4 a b}}\right ) du $$

Now, it should be understood that as $x$ traverses from $0$ to $\infty$, $u$ traverses from $\infty$ down to a min of $2 \sqrt{a b}$ (corresponding to $x \in [0,\sqrt{b/a}]$), then from $2 \sqrt{a b}$ back to $\infty$ (corresponding to $x \in [\sqrt{b/a},\infty)$). Therefore the integral is

$$\frac1{2 a} \int_{\infty}^{2 \sqrt{a b}} du \left (1 - \frac{u}{\sqrt{u^2-4 a b}}\right ) e^{-u} + \frac1{2 a} \int_{2 \sqrt{a b}}^{\infty} du \left (1 + \frac{u}{\sqrt{u^2-4 a b}}\right ) e^{-u} $$

which simplifies to

$$\frac1{a} \int_{2 \sqrt{a b}}^{\infty} du \frac{u}{\sqrt{u^2-4 a b}} e^{-u} = 2 \sqrt{\frac{b}{a}} \int_0^{\infty} dv \cosh{v} \, e^{-2 \sqrt{a b} \cosh{v}}$$

which is then

$$\int_0^{\infty} dx \, e^{-(a x+b/x)} = 2 \sqrt{\frac{b}{a}} K_1\left ( 2 \sqrt{a b}\right )$$

As you say, this is a difficult integral. Provided $\Re(a)>0$ and $\Re(b)>0$, the solution I found is $$\frac{2 \sqrt{b} K_1\left(2 \sqrt{a} \sqrt{b}\right)}{\sqrt{a}}$$ in which appears the modified Bessel function of the second kind.

I obtained the result from a CAS. I have not been able to obtain any analytical form for the antiderivative but, surprizingly, the integral came out !

Of course we need $a,b>0$.

Knowing that the answer is supposed to be what Claude obtained, I can confirm that. By linear change of variables, we may assume $a=1$. Also for convenience take $b = c^2/4$. So now let $$F(c) = \dfrac{1}{c} \int_0^\infty e^{-x - c^2/(4 x)}\ dx$$ I claim $F(c) = K_1(c)$.
A differential equation for $y(c) = K_1(c)$ is $$ y'' + \dfrac{1}{c} y' - \left(1 - \dfrac{1}{c^2}\right) y = 0$$ If we apply that to $F(c)$ and do some simplification, we get $$ F'' + \dfrac{1}{c} F' - \left(1 - \dfrac{1}{c^2}\right) F = \int_0^\infty e^{-x-c^2/(4x)} \dfrac{c^2 - 4 x^2}{4 c x^2}\ dx $$ Using the change of variables $x = c^2/(4t)$, the right side becomes $$ - \int_0^\infty e^{-t-c^2/(4t)} \dfrac{c^2 - 4 t^2}{4 c t^2}\ dt $$ and therefore is $0$. Now the general solution of the differential equation is $A I_1(c) + B K_1(c)$ where $A$ and $B$ are arbitrary constants. $K_1(c)$ is the solution that goes to $0$ as $c \to +\infty$ and is asymptotic to $1/c$ as $c \to 0+$.