Evaluating $\int_0^\infty (\operatorname{E}_n(x)e^x-\frac1{1+x})dx$ 
I want to evaluate $$I_n=\int_0^\infty \left(\operatorname{E}_n(x)e^x-\frac1{1+x}\right)dx=-\psi(n)$$where $\operatorname{E}_n$ denotes exponential integral and $\psi$ denotes polygamma function.  

Attempt
$$I_n=\int_0^\infty\left(\int_1^\infty\frac{e^{xt}e^x}{t^n}dt-\frac1{1+x}\right)dx\\
=\int_0^\infty\left(\int_0^\infty\frac{e^{xt}}{(t+1)^n}dt-\frac1{1+x}\right)dx\\
=\int_0^\infty\left(\int_0^\infty\frac{e^{xt}}{(t+1)^n}-e^{-(x+1)t}\right)dtdx$$
I tried to change the position of two integral signs but failed because of divergency.
 A: Since 
$$\operatorname{E}_n (x) \equiv \int \limits_1^\infty \frac{\mathrm{e}^{-x t}}{t^n} \, \mathrm{d} t \sim \frac{\mathrm{e}^{-x}}{x}\left[1 - \frac{n}{x} +\mathcal{O}\left(\frac{1}{x^2}\right)\right]$$
holds as $x \to \infty$ , $I_n$ exists for $n \in \mathbb{N}$ . We have $\operatorname{E}_n' = - \operatorname{E}_{n-1}$ for $n \in \mathbb{N}$ by definition, so we can integrate by parts to obtain
$$ I_n = \operatorname{E}_n(x) \mathrm{e}^x ~\Bigg \vert_{x=0}^{x=\infty} + \int \limits_0^\infty \left[\operatorname{E}_{n-1}(x) \mathrm{e}^x - \frac{1}{1+x}\right] \, \mathrm{d} x = - \frac{1}{n-1} + I_{n-1} $$
for $n > 1$ . This implies $I_n = I_1 - H_{n-1}$ for $n \in \mathbb{N}$ . In order to evaluate $I_1$ we integrate by parts again, now using the antiderivative $x \mapsto \mathrm{e}^x - 1$ of $x \mapsto \mathrm{e}^x$ , however. We find
$$ I_1 = \operatorname{E}_1(x) (\mathrm{e}^x-1) ~\Bigg \vert_{x=0}^{x=\infty} + \int \limits_0^\infty \left[\frac{\mathrm{e}^{-x}}{x} (\mathrm{e}^x - 1) - \frac{1}{1+x}\right] \, \mathrm{d} x  = \int \limits_0^\infty \left[\frac{1}{1+x} - \mathrm{e}^{-x}\right] \, \frac{\mathrm{d} x}{x} = \gamma$$
from a well-known integral representation of the Euler-Mascheroni constant, whence the desired result follows:
$$ I_n = \gamma - H_{n-1} = - \psi(n) \, , \, n \in \mathbb{N} \, .$$
