This is an example for the asymptotic series. The asymptotic series of the integral is
$$I(\lambda ) \sim \sum_{n=0}^\infty (-1)^n \frac{n!}{\lambda^{n+1}} .$$
It is possible to get a high precision estimate for $I$ when $\lambda $ is moderately large. However, by this approach, we cannot diminish the error to arbitrarily small for given $\lambda$.
By change of variable, $w = e^{-\lambda x }$, we can convert the improper integral to the proper integral
$$ I(\lambda ) = \int_0^1 \frac{d w}{\lambda - \ln w }$$
But there is still some singularity at $w=0$. So how can one calculate the integral to arbitrary precision? Any good reference on numerical integration?
Concerning the antiderivative $$I=\int\frac{e^{-\lambda x}}{x+1}\,dx$$ let $1+x=t$ to make $$I=e^{\lambda} \int \frac{e^{-\lambda t}}{t}\,dt=e^{\lambda}\int\frac{e^{-u}}{u}\,du=e^{\lambda}\,\text{Ei}(-u)=e^{\lambda }\,\text{Ei}(-(x+1) \lambda )$$ So, if $\lambda >0$, then $$\color{blue}{J=\int_0^\infty\frac{e^{-\lambda x}}{x+1}\,dx=-e^{\lambda } \,\text{Ei}(-\lambda )=e^{\lambda }\, \Gamma (0,\lambda )}$$