I have used Mathematica, I must confess, but this integral is known in a closed form
$$
\int_0^\infty dx e^{-kx}[I_0(2x)]^2=\frac{2}{k\pi}K\left(\frac{16}{k^2}\right).
$$
The condition $k>4$ must hold. Here $K(m)$ is an elliptic integral given by
$$
K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-m \sin^2\theta}} = \int_0^1 \frac{dt}{\sqrt{(1-t^2)(1-m t^2)}}.
$$