$\int_0^\infty \frac{\log(1+x)}{x}e^{-\alpha x}dx$ I would like to calculate
$$
F(\alpha)=\int_0^\infty \frac{\log(1+x)}{x}e^{-\alpha x}dx
$$
for $\alpha>0.$ 
Since WolframAlpha provides a complicated expression in terms of hypergeometric functions and the exponential integral $\mathrm{Ei}$ function already for $\alpha=1$, I see little hope for an exact analytical evaluation. Therefore I tried to estimate $F(\alpha)$ at least in some limit, as follows: integrating repeatedly by parts,
$$\begin{aligned}
f(\alpha)\equiv -F'(\alpha)&=\int_0^\infty \log(1+x)e^{-\alpha x}dx\\
&=-\frac{1}{\alpha}\log(1+x)e^{-\alpha x}\bigg|_0^\infty+\frac{1}{\alpha}\int_0^\infty \frac{e^{-\alpha x}}{1+x}dx\\
&=0-\frac{e^{-\alpha x}}{\alpha^2(1+x)}\bigg|_0^\infty-\frac{1}{\alpha^2}\int_0^\infty\frac{e^{-\alpha x}}{(1+x)^2}dx\\
&=\frac{1}{\alpha^2}+\frac{e^{-\alpha x}}{\alpha^3(1+x)^2}\bigg|_0^\infty+\frac{2}{\alpha^3}\int_0^\infty \frac{e^{-\alpha x}}{(1+x)^3}dx\\
&=\frac{1}{\alpha^2}-\frac{1}{\alpha^3}-\frac{2e^{-\alpha x}}{\alpha^4(1+x)^3}\bigg|_0^\infty-\frac{3!}{\alpha^4}\int_0^\infty\frac{e^{-\alpha x}}{(1+x)^4}dx.
\end{aligned}$$
Proceeding by induction,
$$
f(\alpha)=\frac{1}{\alpha^2}\left(1-\frac{1}{\alpha}+\frac{2}{\alpha^2}-\ldots+\frac{(-1)^nn!}{\alpha^n}\right)+(-1)^{n+1}\frac{(n+1)!}{\alpha^{n+2}}\int_0^\infty\frac{e^{-\alpha x}}{(1+x)^{n+2}}dx,
$$
from which
$$\begin{aligned}
\lim_{\alpha\to\infty}\alpha^{n+2}\left[f(\alpha)-\sum_{j=0}^{n-1}\frac{(-1)^jj!}{\alpha^{j+2}}\right]&=(-1)^{n}{n!}+(-1)^{n+1}(n+1)!\lim_{\alpha\to\infty}\int_0^\infty\frac{e^{-\alpha x}}{(1+x)^{n+2}}dx\\
&=(-1)^{n}{n!}
\end{aligned}$$
by dominated convergence. This shows that we have obtained an asymptotic series for $f(\alpha)$ in the limit $\alpha\to\infty$:
$$
f(\alpha)\underset{\alpha\to\infty}{\sim}\sum_{n=0}^\infty \frac{(-1)^nn!}{\alpha^{n+2}}\ .
$$
Even if it is of course divergent, this series is Borel summable: 
$$\boxed{
\sum_{n=0}^\infty \frac{(-1)^nn!}{\alpha^{n+2}}\overset{B}{=}\int_0^\infty dt\, e^{-t}\sum_{n=0}^\infty \frac{(-1)^nn!}{\alpha^{n+2}}\frac{t^n}{n!}\,{\color{red}=}\,\frac{1}{\alpha^2}\int_0^\infty e^{-t(1+1/\alpha)}dt=\frac{1}{\alpha(1+\alpha)}}\ .
$$
It is therefore tempting to state that
$$
f(\alpha)=\int_0^\infty \log(1+x)e^{-\alpha x}dx\overset{?}{=}\frac{1}{\alpha(1+\alpha)}
$$
and hence, integrating in $\alpha$ and setting the integration constant to zero because $F(\alpha)\to0$ as $\alpha\to\infty$ by dominated convergence,
$$
F(\alpha)=\int_0^\infty \frac{\log(1+x)}{x}e^{-\alpha x}dx \overset{?}{=} \log\frac{1+\alpha}{\alpha}.
$$
Strictly speaking, these equalities are wrong, but it appears that they provide reliable numerical approximations to $f(\alpha)$ and $F(\alpha)$ in the $\alpha\to\infty$ limit. 
How do I keep track of the degree of approximation involved in the resummation? This was explicit in the asymptotic series, which is however not as  nice to write down.
EDIT:
Maybe it can be instructive to compare with a similar situation where the exact solution is available. Consider
$$
G(\alpha)=\int_0^\infty\frac{\sin x}{x}e^{-\alpha x}dx
$$
for $\alpha>0$. 
Here
$$
g(\alpha)\equiv-G'(\alpha)=\int_0^\infty \sin x\, e^{-\alpha x}dx=\mathrm{Im}\left[\int_0^\infty e^{-(\alpha-i)x}dx\right]=\mathrm{Im}\left[\frac{1}{\alpha-i}\right]=\frac{1}{1+\alpha^2}
$$
and hence, because $G(\alpha)\to0$ as $\alpha\to\infty$,
$$
G(\alpha)=\frac{\pi}{2}-\arctan \alpha.
$$
Nevertheless we can work out an asymptotic series for $g(\alpha)$ by integrating by parts. The result is
$$
g(\alpha)\underset{\alpha\to\infty}{\sim}\sum_{n=0}^\infty \frac{(-1)^n}{\alpha^{2n+2}}.
$$
In this case however the series converges to the exact result $g(\alpha)=1/(1+\alpha^2)$ for $\alpha>1$.
EDIT 2: as was pointed out by Sangchul Lee in the comments, the derivation in the box contains a mistake
$$
\sum_{n=0}^{^\infty}\frac{(-1)^n n!}{\alpha^{n+2}}\overset{B}{=}\int_0^\infty \frac{e^{-t}}{\alpha(t+\alpha)}dt\,{\color{red}\neq}\,\frac{1}{\alpha(1+\alpha)},
$$
which is one of the integral representations of $f(\alpha)$ already appearing in the second equation. Therefore, Borel resummation gives the exact result. What I still don't know is how the wrong expressions I derived give a reasonable approximation for large $\alpha$.
 A: I'll provide a somehow different approach to approximate the integral 
$$F(\alpha) = \int_0^\infty \frac{\log(1+x)}{x}e^{-\alpha x}dx$$
Using integration by parts 
$$F(\alpha) = \int_0^\infty \frac{\log(1+x)}{x}e^{-\alpha x}dx=-\alpha\int^\infty_0 e^{-\alpha x}\mathrm{Li}_2(-x)\,dx$$
For the record this integral appears in Lewis book
$$\int^\infty_0 e^{-\alpha x} \mathrm {Li}_2 (-x) \, dx = \frac {1}{\alpha}\int^\infty_{\alpha}\frac {e^x}{x}\mathrm {Ei}(-x)\, dx$$
Hence we have 
$$F(\alpha) = \int^\infty_{\alpha}\frac {e^x}{x}\mathrm {E}_1(x)\, dx$$
Now use the approximation 
$$\frac{1}{2}e^{-x}\log\left( 1+\frac{2}{x}\right)<\mathrm{E}_1(x) < e^{-x}\log\left( 1+\frac{1}{x}\right)$$
Wiki picture showing the tightness of this bound 

Hence we have 
$$\frac{1}{2} \int^\infty_{\alpha}\frac{\log\left( 1+\frac{2}{x}\right)}{x}\,dx<F(\alpha) <  \int^\infty_{\alpha}\frac{\log\left( 1+\frac{1}{x}\right)}{x}\,dx$$
This can be rewritten as 

$$ -\frac{1}{2}\mathrm{Li}_2\left( -\frac{2}{\alpha}\right)< F(\alpha)
 < -\mathrm{Li}_2\left( -\frac{1}{\alpha}\right)$$

A plot for 30 points for $\alpha$ showing the upper and lower bounds

A plot for 10 points for $\alpha$ showing the upper and lower bounds

A scatter plot for 10 points 

A: Behavior of the Laplace transform $\mathcal{L}\{g\}(s)$ for large parameter $s$ is intimately related to the near-zero behavior of the function $g$ being transformed. In particular, the leading order of $\mathcal{L}\{g\}(s)$ as $s \to \infty$ is very robust and roughly depends only on the value $g(0)$. That is a reason why your computation still gives a good approximation.

Next, let me derive an asymptotic expansion for $F(\alpha)$. We begin by splitting the integral into two parts
$$F(\alpha)
= \int_{0}^{\epsilon} \frac{\log(1+x)}{x} e^{-\alpha x} \, dx + \int_{\epsilon}^{\infty} \frac{\log(1+x)}{x} e^{-\alpha x} \, dx. $$
Applying the Cauchy-Schwarz inequality to the latter term, we have
\begin{align*}
\left|\int_{\epsilon}^{\infty} \frac{\log(1+x)}{x} e^{-\alpha x} \, dx\right|
&\leq \bigg( \int_{\epsilon}^{\infty} \frac{\log^2(1+x)}{x^2}\, dx \bigg)^{1/2}\bigg( \int_{\epsilon}^{\infty} e^{-2\alpha x} \, dx \bigg)^{1/2} \\
&\leq C\alpha^{-1/2}e^{-\epsilon\alpha}
\end{align*}
and thus the latter term only contributes to an exponentially decaying error. On the other hand, if $\epsilon < 1$ and $N$ is a positive integer, then
\begin{align*}
\int_{0}^{\epsilon} \frac{\log(1+x)}{x} e^{-\alpha x} \, dx
&= \sum_{n=1}^{N} \frac{(-1)^{n-1}}{n} \int_{0}^{\epsilon} x^{n-1} e^{-\alpha x} \, dx \\
&\qquad + \int_{0}^{\epsilon} \bigg( \frac{\log(1+x)}{x} - \sum_{n=1}^{N} \frac{(-1)^{n-1}}{n}x^{n-1} \bigg) e^{-\alpha x} \, dx \\
&= \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n} \frac{1}{\alpha^n} \int_{0}^{\alpha \epsilon} x^{n-1} e^{-x} \, dx \\
&\qquad + \mathcal{O}\bigg( \int_{0}^{\infty} x^N e^{-\alpha x} \, dx \bigg).
\end{align*}
Here, for each fixed $n$, it is not hard to check that
$$\int_{0}^{\alpha \epsilon} x^{n-1} e^{-x} \, dx = (n-1)! + \mathcal{O}\big( (\epsilon\alpha)^{n-1}e^{-\epsilon\alpha} \big). $$
Putting altogether, for each fixed $N$ we have an asymptotic expansion

$$ F(\alpha) = \sum_{n=0}^{N-1} \frac{(-1)^n n!}{n+1} \frac{1}{\alpha^{n+1}} + \mathcal{O}\bigg(\frac{1}{\alpha^{N+1}}\bigg) \qquad \text{as } \alpha \to \infty $$

It is not surprising that the first term is also obtained formally by integrating OP's asymptotic expansion for $f(\alpha)$ as well.

For the near-zero behavior of $F$, we begin from the expression @Zaid Alyafeai derived:
$$ F(\alpha)
= -\alpha \int_{0}^{\infty} \operatorname{Li}_2(-x) e^{-\alpha x} \, dx
= -\int_{0}^{\infty} \operatorname{Li}_2(-x/\alpha) e^{-x} \, dx. $$
Utilizing the identity
$$ -\operatorname{Li}_2(-z) = \zeta(2) + \frac{1}{2}\log^2 z + \operatorname{Li}_2 (-1/z), $$
we find that
\begin{align*}
F(\alpha)
&= \int_{0}^{\infty} \left( \zeta(2) + \frac{1}{2}\log^2(x/\alpha) + \operatorname{Li}_2(-\alpha/x) \right) e^{-x} \, dx \\
&= \frac{1}{2}\log^2\alpha + \gamma \log\alpha + \frac{1}{2}(\gamma^2 + 3\zeta(2)) + \int_{0}^{\infty} \operatorname{Li}_2(-\alpha/x) e^{-x} \, dx.
\end{align*}
We claim that the last integral vanishes as $\alpha \to 0$. Indeed, from the estimate $|\operatorname{Li}_2(-x^{-1})| \sim x^{-1}$ as $x \to \infty$, it is possible to prove that
$$ \int_{0}^{\infty} \operatorname{Li}_2(-\alpha/x) e^{-x} \, dx
= \alpha \int_{0}^{\infty} \operatorname{Li}_2(-1/x) e^{-\alpha x} \, dx
= \mathcal{O}\big(\alpha \log (1/\alpha) \big). $$
Therefore

$$F(\alpha) = \frac{1}{2}\log^2\alpha + \gamma \log\alpha + \frac{1}{2}(\gamma^2 + 3\zeta(2)) + \mathcal{O}\big(\alpha \log (1/\alpha) \big) \qquad \text{as } \alpha \to 0^+. $$

In fact a more detailed analysis on the error term is available. Let
$$ g(\alpha) = \int_{0}^{\infty} \operatorname{Li}_2(-\alpha/x) e^{-x} \, dx. $$
Then $g$ defines a continuous function on $[0,\infty)$ which is differentiable on $(0,\infty)$. Differentiating under the integral sign,
$$ g'(\alpha)
= - \int_{0}^{\infty} \frac{1-e^{-x}}{x(x+\alpha)}\, dx
= - \frac{1}{\alpha} \int_{0}^{\infty} \left( \frac{1}{x} - \frac{1}{x+\alpha} \right)(1-e^{-x}) \, dx. $$
The last integral exhibits cancellation of singularity at infinity. In order to analyze this effect, we work with the truncated integral.
\begin{align*}
&\int_{0}^{R} \left( \frac{1}{x} - \frac{1}{x+\alpha} \right)(1-e^{-x}) \, dx \\
&= \int_{0}^{\alpha} \frac{1-e^{-x}}{x} \, dx + \int_{\alpha}^R \frac{1-e^{-x}}{x} \, dx - \int_{\alpha}^{R+\alpha} \frac{1-e^{-(x-\alpha)}}{x} \, dx \\
&= \int_{0}^{\alpha} \frac{1-e^{-x}}{x} \, dx + \log \left(\frac{R}{R+\alpha}\right) - \int_{\alpha}^R \frac{e^{-x}}{x} \, dx + e^{\alpha}\int_{\alpha}^{R+\alpha} \frac{e^{-x}}{x} \, dx.
\end{align*}
Taking limit as $R\to\infty$, we find that
$$ g'(\alpha) = - \frac{1}{\alpha} \int_{0}^{\alpha} \frac{1-e^{-x}}{x} \, dx -
 \frac{e^{\alpha} - 1}{\alpha} E_1(\alpha), \tag{1} $$
where $E_1(\alpha) = \int_{\alpha}^{\infty} e^{-x}/x \, dx$ is a variant of the exponential integral. It is well-known that $E_1(\alpha)$ has the following expansion
$$ E_1(\alpha) = -\gamma - \log \alpha - \sum_{n=1}^{\infty} \frac{(-1)^n}{n!n} \alpha^n $$
which is valid for all $\alpha \in \Bbb{C} \setminus (-\infty, 0]$. Plugging this back, $g'(\alpha)$ take the form
$$ g'(\alpha) =  \frac{e^{\alpha} - 1}{\alpha} \log \alpha + \text{[entire function in $\alpha$]}. $$
Mathematical tells that this entire-function part has a neat series expansion, yielding
$$ g'(\alpha) = \frac{e^{\alpha} - 1}{\alpha} \log \alpha - \sum_{n=1}^{\infty} \frac{-\gamma + H_n}{n!}x^{n-1}. $$
I haven't checked this by myself, but I expect that this is not hard to verify by expanding everything in $\text{(1)}$. Integrating, end up with the following series expansion

$$F(\alpha) = \frac{1}{2}\log^2\alpha + \gamma \log\alpha + \frac{1}{2}(\gamma^2 + 3\zeta(2)) + \sum_{n=1}^{\infty} \frac{\alpha^n}{n!n} \left(\gamma + \log \alpha - H_n - \frac{1}{n} \right). $$

For instance, the following is a numerical computation of $F(2)$ using both numerical integration and the formula above:


I suspect that more systematic approach is available for both directions, but I am not sure how to proceed.
A: Elaborating slightly on the answer by Sangchul Lee,  I would like to elucidate why my (wrong) computation still gives a reasonable approximation.
The mistake was made in writing
$$
f(\alpha)=\int_0^\infty \frac{e^{-\alpha t}}{\alpha(t+1)}dt\,{\color{\red}=}\,\frac{1}{\alpha(1+\alpha)}.
$$
Since $e^{-\alpha t}$ attains its global maximum on the integration region when $t=0$, this integral may be approximated for large $\alpha$ by Laplace's method:
$$
f_\epsilon(\alpha) = \frac{1}{\alpha}\sum_{n=0}^\infty (-1)^n\int_0^{\epsilon}t^n e^{-\alpha t}dt;
$$
replacing $\epsilon$ with $\infty$ now only introduces exponentially small errors so that
$$
f(\alpha)\underset{\alpha\to\infty}\sim \sum_{n=0}^\infty \frac{(-1)^nn!}{\alpha^{n+2}}=\frac{1}{\alpha^2}-\frac{1}{\alpha^3}+\mathcal O(\alpha^{-4}).
$$
By chance, for large $\alpha$,
$$
\frac{1}{\alpha(1+\alpha)}=\frac{1}{\alpha^2}\frac{1}{1+1/\alpha}=\frac{1}{\alpha^2}\left(1-\frac{1}{\alpha}+\mathcal O(\alpha^{-2})\right)=\frac{1}{\alpha^2}-\frac{1}{\alpha^3}+\mathcal O(\alpha^{-4}).
$$
So the mistake done in the computation gives the right answer up to order $\alpha^{-4}$ because it is the sum of a geometric series which agrees with the right asymptotic expansion for $f(\alpha)$ precisely up to that order.
