Integral with Bessel function and hypergeometric function ${}_2F_2$: explicit expression for these polynomials? This question follows this one, where the general problem has apparently no simpler form than the integral one. I focus now on the limit case:
\begin{align}
\int_0^T e^{-x}\frac{nI_n(x)}{x}dx=\int_0^T e^{-x}\frac{I_{n-1}(x)-I_{n+1}(x)}{2}dx
\end{align}
where $I_n(x)$ is the modified Bessel function of the first kind, and $n$ is an integer.
I can get the following result from the series expansion of the modified Bessel function:
\begin{align}
\int_0^T e^{-x}\frac{I_{n}(x)}{2}dx&=\frac{(T/2)^{n+1}}{(n+1)!}{}_2F_2[\{n+\frac{1}{2},n+1\};\{2n+1,n+2\};-2T] \\
\int_0^T e^{-x}\frac{nI_n(x)}{x}dx&=\frac{(T/2)^{n}}{n!}{}_2F_2[\{n+\frac{1}{2},n\};\{2n+1,n+1\};-2T]
\end{align}
I have some trouble to transforming this expression. I know from the answer to this post that there exist some polynomials of order $n$, $P_n(T)$ and $Q_n(T)$, such that:
\begin{align}
\int_0^T e^{-x}\frac{nI_n(x)}{x}dx=\Big(1+\frac{e^{-T}}{T^{n-1}}\big(P_{n-1}(T)I_0(T)+Q_{n-1}(T)I_1(T)\big)\Big)
\end{align}
which looks somehow to the result of an integration by parts. However I don't manage to use the derivation relation for the Bessel function, or the identities (especially number 2) with the hypergeometric function here to obtain the expression of the polynomials.
Furthermore, I have some hints (from numerical evaluations of the formula of a related model) that this could be written in terms of error functions, powers and exponentials. I have troubles to use the previous identities.
What would be an effective method to get more informations from this integral?
 A: From DLMF, one has, for $n=0,1,2...$
$$\int_{0}^{x}e^{-t}I_{n}\left(t\right)\mathrm{d}t=xe^{-x}(I_{0}\left(x\right)+I%
_{1}\left(x\right))+n(e^{-x}I_{0}\left(x\right)-1)+2e^{-x}\sum_{k=1}^{n-1}(n-k%
)I_{k}\left(x\right)$$
which may be obtained by recurrence, I guess. Then
\begin{align}
F_n&=\frac{1}{2}\int_0^T e^{-x}[I_{n-1}(x)-I_{n+1}(x)]\,dx\\
&=1-e^{-T}\left[I_{0}\left( T \right)+I_{n}\left( T \right)+2\sum_{k=1}^{n-1}I_{k}\left(T\right)\right]  \label{eq:1}\tag{1}
\end{align} 
This expression can be checked by taking the derivative of the rhs, using the recursion relation for the Bessel function $I'_n(T)=\frac{1}{2}\left(I_{n-1}(T)+I_{n+1}(T)\right)$:
$$F'_n=\frac{1}{2}e^{-T}[I_{n-1}(T)-I_{n+1}(T)]$$
Now, using the recurrence relation for the Bessel functions
\begin{equation}
I_{k}\left( T \right)=I_{k-2}\left( T \right)-\frac{2(k-1)}{T}I_{k-1}\left( T \right)
\end{equation} 
one deduce that
\begin{equation}
I_k\left( T \right)=T^{1-k}\left[A_{k-1}\left( T \right)I_0\left( T \right)+B_{k-1}\left( T \right)I_1\left( T \right)\right]
\end{equation} 
where $A_{k-1}$ and $B_{k-1}$ are polynomials of degree $k-1$.
Finally, as expected
\begin{equation}
F_n=1+\frac{e^{-T}}{T^{n-1}}\left[P_{n-1}\left( T \right)I_0\left( T \right)+Q_{n-1}\left( T \right)I_1\left( T \right)\right]
\end{equation} 
Eq. (\ref{eq:1}) can be evaluated by Maple using ``simplify'' command giving the polynomials obtained in the previous post Closed form for $\int_0^T e^{-x}\frac{I_n(\alpha x)}{x}dx$
Edit (22/08/'17): The polynomials can be expressed using analogous of the Lommel polynomials for the modified Bessel functions. In Lemma 5 of this paper,
$$
I_k(T)=\left( \frac{2}{T} \right)^k r_k^0\left( \frac{T^2}{4} \right)I_0(T)+\left( \frac{2}{T} \right)^{k-1} s_k^0\left( \frac{T^2}{4} \right)I_1(T)
$$
Expressions for the polynomials $r_k^0(x)$ and $s_k^0(x)$ are given either in terms of the Lommel polynomials or explicitly.
A: In addition to this excellent answer, I would like to answer this secondary question:

Furthermore, I have some hints (from numerical evaluations of the formula of a related model) that this could be written in terms of error functions, powers and exponentials.

I have the answer for this. The study of this formula was an asymptotic limit of the expression of this other related question:
\begin{align}
\int_0^T e^{-x}\frac{nI_n(x)}{x}dx=\lim_{\alpha \rightarrow1}\int_0^T e^{-x}\frac{nI_n(\alpha x)}{x}dx
\end{align}
and $\alpha \rightarrow1$ in my specific model also meant (not explicitly here) that $x$ and $n$ are large. Then I needed to take the following asymptotic formula using Hankel's expansion of the Bessel function:
\begin{align}
e^{-x}I_n(x)&=\frac{1}{\sqrt{2\pi x}}\sum_{k=0}^\infty (-1)^k \frac{\prod_{l=1}^k \big((2n)^2-(2l-1)^2\big)}{k!\,(8x)^k} \sim \frac{e^{-\frac{n^2}{2x}}}{\sqrt{2\pi}}
\end{align}
The gaussian function appears and the error function too once integrated. This is a continuous limit for the discrete integers $n$.
