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$.