Question: $$I=\large \int\_{0}^{\infty}\left(x^2-3x+1\right)e^{-x}\ln^3(x)\mathrm dx$$

$$e^{-x}=\sum\_{n=0}^{\infty}\frac{(-x)^n}{n!}$$

$$I=\large \sum\_{n=0}^{\infty}\frac{(-1)^n}{n!}\int\_{0}^{\infty}\left(x^2-3x+1\right)x^n\ln^3(x)\mathrm dx$$

$$J=\int \left(x^2-3x+1\right)x^n\ln^3(x)\mathrm dx$$

We can evaluate $J$ by integration by parts but problem, the limits does not work.

**How to evaluate integral $I?$**

Answer: 
It is pretty straightforward to differentiate three times
$$ \int\_{0}^{+\infty}(x^2-3x+1)x^s e^{-x}\,dx = s^2\,\Gamma(s+1)$$
then consider the limit as $s\to 0^+$. The final outcome is $\color{red}{-6\,\gamma}=6\,\Gamma'(1)$ since $s^2\,\Gamma(s+1)$ clearly has a zero of order $2$ at the origin.
The form of the integral takes on the form of the product of a power function, $e^{-x}$, and an integer power of a logarithm. This allows us to consider the integral

$$ \int\_{0}^{\infty}x^{s-1+\epsilon}\,e^{-x}\,\mathrm{d}x = \int\_{0}^{\infty}x^{s-1}e^{\epsilon\ln x}\,e^{-x}\,\mathrm{d}x = \sum\_{n=0}^{\infty}\frac{\epsilon^{n}}{n!}\int\_{0}^{\infty}x^{s-1}\,e^{-x}\,\ln^{n}x\,\mathrm{d}x. $$

Using the gamma function, this integral is also

$$ \int\_{0}^{\infty}x^{s-1+\epsilon}\,e^{-x}\,\mathrm{d}x = \Gamma(s + \epsilon). $$

Using the recursion relation, we may write this in terms of $\Gamma(1+\epsilon)$, which we have a Taylor series expansion for. For example, with $x^{2}$ in the integrand, $s = 3$, so we get $\Gamma(3+\epsilon) = (2+\epsilon)(1+\epsilon)\Gamma(1+\epsilon)$.

We have rewritten our integral as a coefficient in a series expansion. Our integral needs the coefficient of $\epsilon^{3}$ to get the integral with $\ln^{3}x$. We therefore have to obtain the corresponding $\epsilon^{3}$ coefficient in the expansion of the gamma function $\Gamma(1+\epsilon)$.

The usual way of tackling this integral would be to consider the three terms separately and then keep the appropriate terms, but we can do better by evaluating the gamma function terms directly first, which contribute in the following manner:

$$\begin{align} x^{2} - 3x + 1 &\to \Gamma(3 + \epsilon) - 3\Gamma(2+\epsilon) + \Gamma(1+\epsilon) \\
&\to \left((2+\epsilon)(1+\epsilon) - 3(1+\epsilon) + 1\right)\Gamma(1+\epsilon) \\
&\to \epsilon^{2}\,\Gamma(1 + \epsilon). \end{align}$$

Since all that is left is an $\epsilon^{2}$ out in front, we only need keep terms up to first order in $\epsilon$ in the expansion of the gamma function (where $\gamma$ is the Euler-Mascheroni constant and $\zeta(s)$ is the Riemann zeta function)

$$\begin{align} \ln\Gamma(1 + \epsilon) &= -\gamma\epsilon + \sum\_{k=2}^{\infty}\frac{(-1)^{k}\zeta(k)}{k}\epsilon^{k} \\
\Gamma(1 + \epsilon) &\approx e^{-\gamma\epsilon} \approx 1 - \gamma\epsilon. \end{align}$$

Our coefficient is therefore $-\gamma$. Remembering to multiply by $3! = 6$ to account for the factorial in the original series expansion, our answer is

$$ \int\_{0}^{\infty}\left(x^{2} - 3x + 1\right)e^{-x}\ln^{3}x\,\mathrm{d}x = -6\gamma. $$