I'm trying to solve the definite integral
$$I_1=\int_0^\infty\exp \left( -\sqrt{x^2+y^2}\, \right)\left(\frac{1}{x^2+y^2}+\frac{1}{\sqrt{x^2+y^2}}\right)dx,$$
with $y>0$ and is obviously symmetric (so boundaries can be changed to $\pm\infty$).
I have some success solving parts of the integral
$$I_2=\int_0^\infty\exp\left(-\sqrt{x^2+y^2}\,\right)\frac{1}{\sqrt{x^2+y^2}} \, dx = \int_y^\infty \exp(-r)\frac{1}{\sqrt{r^2-y^2}} \, dr=K_0(y),$$
where $K_0$ is the modified Bessel Funktion of the second kind (or MacDonald function). I'm not sure if it's easier solve the remaining part or the whole integral at once with some simplification, that's why I'm giving to full integral above. I tried integration by parts and the residue theorem for $I_1$ and $I_1-I_2$, but never arrived at something useful.
From literature I know as well
$$ \begin{align} I_3 & =\int_0^\infty y\exp(-\sqrt{x^2+y^2})\left(\frac{1}{x^2+y^2}+\frac{1}{(x^2+y^2)^{3/2}}\right)dx \\[8pt] & = \int_0^\infty \frac{1}{y}\exp(-\sqrt{x^2+y^2})\,dx=K_1(y), \end{align} $$
which I honestly don't understand fully. It looks like the Integral is somehow rewritten by integration by parts, but trying that I never arrive at something similar or useful. I expect there is a trick due to the boundaries and symmetry of the integral...
In summary I need help solving $I_1$ or $I_1-I_2$, or explaining $I_3$.
Thank you very much.
Funnily I found the answer shortly after I posted this question... I have been working on that integral on/off for two weeks.
The result is
$$ I_1 = \frac{\pi}{2} \left(\frac{1}{y} - K_0(y)L_{-1}(y) - K_1(y)L_0(y)\right) + K_0(y) \; . $$
A coworker and I got the result by rewriting the integral, taking the derivative of the integral and solving the resulting differential equation. If there is any interest I can describe it in more detail...