Evaluating $I(x)=\int_{0}^{\infty}\frac{e^{-xy}}{y^2+a^2}dy$ I examine currently this integral:   
$$I(x)=\int_{0}^{\infty}\frac{e^{-xy}}{y^2+a^2}dy;x\geqslant0$$ where $x$ and $a$ are real.  
It seems that the integral has no closed form in terms of elementary functions. But perhaps it has a closed form in terms of special functions?  
The end result should be series expansion of $I(x)$ at $x=0$  
Thanks!
 A: Note that $$I''(x)=\int_0^\infty \frac{y^2 e^{-xy}}{y^2+a^2}dy$$
so that $$\tag{1} I''(x)+a^2 I(x)= \int_0^\infty e^{-xy} dy=\frac 1x$$
Solutions of the homogenous equation are $I(x)=\sin(ax)$ and $I(x)=\cos(ax)$ but specific solutions will require (as noticed by Norbert) sine and cosine integrals as we will show using variation of constants :
$$\tag{2} I(x)=c\;\sin(ax)$$
$$I'(x)=c'\sin(ax)+c\;a\cos(ax)$$
$$I''(x)=c''\sin(ax)+2\,c'a\cos(ax)-c\;a^2\sin(ax)$$
$$I''(x)+a^2 I(x)=\frac 1x=c''\sin(ax)+2\,c'a\cos(ax)$$
Let's set $\;b:=c'$ then we want $$b'\sin(ax)+2\,b\,a\cos(ax)=\frac 1x$$
multiply this by $\,\sin(ax)$ to get :
 $$b'\sin(ax)^2+b\,2\,a\sin(ax)\cos(ax)=\frac {\sin(ax)}x$$
 $$b'\sin(ax)^2+b \frac d{dx} \sin(ax)^2=\frac {\sin(ax)}x$$
 $$\frac d{dx} (b\;\sin(ax)^2)=\frac {\sin(ax)}x$$ 
or $$c'=b=\frac 1{\sin(ax)^2} \left(C_0+\int \frac {\sin(ax)}{ax} d(ax)\right)=\frac {C_0+\rm{Si(ax)}}{\sin(ax)^2}$$
and $c$ will be : 
$$c=C_1+\int \frac {C_0+\rm{Si(ax)}}{\sin(ax)^2} dx$$
but (int. by parts and since $\cot(x)'=-\frac 1{\sin(x)^2}$, $\rm{Si}'(x)=\frac {\sin(x)}x$ and $\rm{Ci}'(x)=\frac {\cos(x)}x$) : 
$$\int \frac {\rm{Si}(ax)}{\sin(ax)^2} dx=\frac 1a\left[-\rm{Si}(ax)\cot(ax)\right]+\int \frac {\sin(ax)}{x} \cot(ax)dx$$
$$\int \frac {\rm{Si}(ax)}{\sin(ax)^2} dx=-\frac 1a\left[\rm{Si}(ax)\cot(ax)\right]+\int \frac {\cos(ax)}{x}dx$$
so that (for $C_0=C_1=0$) $c$ will simply be :
$$\tag{3} c=\frac{\rm{Ci}(ax)-\rm{Si}(ax)\cot(ax)}a$$
multiplying $c$ by $\,\sin(ax)$ in $(2)$ we get the specific solution :
$$\tag{4} I(x)=\frac 1a\left(\rm{Ci}(ax)\sin(ax)-\rm{Si}(ax)\cos(ax)\right)$$
with the general solution of the O.D.E. given by :
$$\tag{5} I(x)=C\,\sin(ax)+D\,\cos(ax)+\frac 1a\left(\rm{Ci}(ax)\sin(ax)-\rm{Si}(ax)\cos(ax)\right)$$
You may use the specific case $I(0)=\left[\frac {\arctan(ax)}a\right]_0^\infty=\frac {\pi}
{2a}$ and the derivative $I'(0)=-\log(a)$ to find Norbert's expression :
$$\tag{6}\boxed{\displaystyle I(x)=\frac {\pi\cos(ax)}{2a}+\frac 1a\left(\rm{Ci}(ax)\sin(ax)-\rm{Si}(ax)\cos(ax)\right)}$$
For numerical evaluation perhaps that the DLMF link or A&S' book will be helpful.
I got following series expansion of $\displaystyle (a\;I(x))$ (as explained by J.M. there is a logarithmic term coming from the $\rm{Ci}$ function that can't be expanded at $0$) :
$$\left[(\gamma+\ln(ax))\sin(ax)+\frac {\pi}2-(ax)-\frac{\pi}4 (ax)^2+\frac {11}{36}(ax)^3+\frac{\pi}{48}(ax)^4-\frac{137}{7200}(ax)^5+\rm{O}\left((ax)^6\right)\right]$$
(with $\gamma\approx 0.5772156649$ the Euler constant)
To get more terms you may use following expansions :
$$\rm{Ci}(z)=\gamma+ \ln(z) +\sum_{n=1}^\infty \frac{(-1)^nz^{2n}}{(2n)(2n)!}$$
$$\rm{Si}(z)=\sum_{n=0}^\infty \frac{(-1)^nz^{2n+1}}{(2n+1)(2n+1)!}$$
P.S. It is quite interesting to notice that there is a reference to nearly the same integral in this other recent S.E. thread, the paper is Coffey 2012 'Certain logarithmic integrals, including solution of Monthly problem #tbd, zeta values, and expressions for the Stieltjes constants' where equation $(4.4)$ reads $\ \displaystyle I(k)\equiv \int_0^\infty \frac{e^{-(k-1)v}}{v^2+\pi^2}\,dv$.
