Parece que la integral se puede expresar mediante McRobert del E-función, que es una función hipergeométrica generalizada. Deje que la función definida como
\begin{equation}
I(X)=\int_0^X e^{\frac{1}{x}+\frac{1}{X-x}} x^{b-a-1}(X-x)^{a-1}\,dx \end{equation}
su valor en $X=1$ se encuentra. Puede ser considerado como una convolución, su transformada de Laplace es así
\begin{equation}
\mathcal{L}\left[I;X\to p\right]=F_{b-a-1}(p)F_{a-1}(p)
\end{equation}
donde
\begin{equation}
F_\alpha(p)=\int_0^\infty x^\alpha e^{-px-\tfrac{1}{x}}\,dx
\end{equation}
De Ederlyi TI, p.146 (4.5.29),
\begin{equation}
F_\alpha(p)=2p^{-\frac{\alpha+1}{2}}K_{\alpha+1}\left( 2\sqrt{p} \right)
\end{equation}
$K_\alpha(.)$ es una función Bessel modificada. Entonces
\begin{equation}
\mathcal{L}\left[I;X\to p\right]=4p^{-\frac{b}{2}}K_{b-a}\left( 2\sqrt{p} \right)K_{a}\left( 2\sqrt{p} \right)
\end{equation}
Ahora, usando la transformada inversa de Laplace del producto de dos funciones de Bessel modificadas derivados de este documento,
\begin{equation}
I(1)=\frac{2^{b-2a}}{\sqrt{\pi}}\sum_{i,-i}\frac{1}{i}\left[E\left( \left.
\begin{array}{l}
b+1,a+1,b-a+1,1\\
\tfrac{b}{2}+1,\tfrac{b+3}{2}
\end{array}
\right|4e^{i\pi} \right)\right]
\end{equation}
donde la suma se toma sobre y debajo de la rama cortada $\Re z<0$ de la McRobert E-función $E\left( \left.
\begin{array}{l}
u_p\\
v_q
\end{array}
\right|z \right)$. This expression may probably be simplified, at least for some specific values of $$ and $b$ e se expresa en términos de la clásica hipergeométrica generalizada funciones o en términos de Meijer de la función G.
Comentario 1: en el Ederlyi tabla (61) p.285, de un caso particular se puede obtener:
$$ \mathcal{L}^{-1}\left[\sqrt{p}K_{\nu+\tfrac{1}{2}}(2\sqrt{p})K_{\nu-\tfrac{1}{2}}(2\sqrt{p})\right]=\frac{1}{4}\sqrt{\frac{\pi}{2}}e^{-2}W_{\tfrac{1}{2},\nu}(4)$$
donde $W_{\mu,\nu}(.)$ es la función de Whittaker. La elección de $\nu=-1/2,b=-1,a=-1$, se trata de
$$I(-1,-1)=\int_0^1 e^{\frac{1}{x}+\frac{1}{1-x}} x^{-1}(1-x)^{-2}\,dx=\sqrt{\frac{\pi}{2}}e^{-2}W_{\tfrac{1}{2},-\tfrac{1}{2}}(4) $$
Sin embargo, numéricamente este resultado es exacto si la r.h.s. se multiplica por un factor de $\sqrt 2$. Errata en la tabla o error en mi derivación?
Observación 2: el Cambio de $x\to 1-x$ en la integral da la identidad
$$I(a,b)=I(b-a,b)$$
Por otra parte, por la relación $(1-x)^a=(1-x)^{a-1}-x(1-x)^{a-1}$, tenemos la propiedad
$$I(a+1,b)=I(a,b-1)-I(a,b)$$
Esta expresión puede ser generalizado, por el uso de la expansión binomial de $(1-x)^{a+k}$. Con estos recurrencia de la relación y los casos particulares de @Nemo y a partir de esta respuesta, podemos acceder a muchos de los casos especiales de la integral.