Evaluate: $\int_{0}^{\infty}\left(x^2-3x+1\right)e^{-x}\ln^3(x) dx$ $$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?$
 A: Here's a way that only requires the identity $\Gamma'(1) = - \gamma$ ,  since all derivatives of higher order cancel:
\begin{align}
I &=\int \limits_0^\infty (x^2 - 3x+1) \ln^3 (x) \mathrm{e}^{-x} \, \mathrm{d} x \\
&= \left[\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 3 \frac{\mathrm{d}}{\mathrm{d}t}+1 \right] \int \limits_0^\infty \ln^3 (x) \mathrm{e}^{- t x} \, \mathrm{d} x ~\Bigg\vert_{t=1}\\
&= \left[\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 3 \frac{\mathrm{d}}{\mathrm{d}t}+1 \right] \frac{1}{t} \int \limits_0^\infty \left[\ln^3 (y) - 3 \ln^2(y) \ln(t) + 3 \ln(y) \ln^2(t) - \ln^3 (t)\right] \mathrm{e}^{-y} \, \mathrm{d} y ~\Bigg\vert_{t=1}\\ 
&= \left[\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 3 \frac{\mathrm{d}}{\mathrm{d}t}+1 \right] \frac{1}{t}  \left[\Gamma'''(1) - 3 \Gamma''(1) \ln(t) + 3 \Gamma'(1) \ln^2(t) - \ln^3 (t)\right] ~\Bigg\vert_{t=1} \\
&= 2 \Gamma'''(1) + 6 \Gamma''(1) + 3 \Gamma''(1) + 6 \Gamma'(1) - 3 \Gamma'''(1) - 9 \Gamma''(1) + \Gamma'''(1) \\
&= 6 \Gamma'(1)\\
&= - 6 \gamma \, .
\end{align}
In fact, integration by parts yields the following generalisation:
\begin{align}
\gamma &=  \int \limits_0^\infty (-\ln (x))  \mathrm{e}^{-x} \, \mathrm{d} x 
 = \int \limits_0^\infty \frac{-\ln (x)}{x} x \mathrm{e}^{-x} \, \mathrm{d} x \\
&= \int \limits_0^\infty (-\ln (x))^2 \frac{1-x}{2} \mathrm{e}^{-x} \, \mathrm{d} x \\
&= \int \limits_0^\infty (-\ln (x))^3 \frac{x^2 - 3x +1}{6} \mathrm{e}^{-x} \, \mathrm{d} x \\
&= \dots \, \\
&= \int \limits_0^\infty (-\ln (x))^{n+1} \frac{p_n (x)}{(n+1)!} \mathrm{e}^{-x} \, \mathrm{d} x \, .
\end{align}
The polynomials $p_n$ are defined recursively by $p_0(x) = 1$ and 
$$p_n (x) = \mathrm{e}^{x} \frac{\mathrm{d}}{\mathrm{d}x} \left(x p_{n-1} (x) \mathrm{e}^{-x}\right) \, , \, n \in \mathbb{N} \, ,$$
for $x \in \mathbb{R}$ . The exponential generating function
$$ \sum \limits_{n=0}^\infty \frac{p_n(x)}{n!} t^n = \mathrm{e}^{t+x(1-\mathrm{e}^t)}$$
can actually be computed from a PDE and it turns out that the polynomials are given by
$$p_n(x) = \frac{B_{n+1}(-x)}{-x} \, , \, x \in \mathbb{R} \, , \, n \in \mathbb{N}_0 \, , $$
where $(B_k)_{k \in \mathbb{N}_0}$ are the Bell or Touchard polynomials.
A: 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.
A: Though it's not good practice, we will ignore any constants when evaluating indefinite integrals, as the end result is definite.
Step I:
Let's start with the integral $$G_1=\int\ln x\,dx=x(\ln x-1).$$ Integration by parts gives $$G_2=\int\ln^2x\,dx=[x\ln x(\ln x-1)]-\int(\ln x-1)\,dx=x(\ln^2x-2\ln x+2)$$ and similarly, $$G_3=\int\ln^3x\,dx=x(\ln^3x-3\ln^2x+6\ln x-6)$$ 
Step II:
Consider the integral $$H_1=\int\frac{\ln x}{e^x}\,dx.$$ Integration by parts gives $$H_1=\left[\frac{G_1}{e^x}\right]+\int\frac{x(\ln x-1)}{e^x}\,dx\implies I_1=\int\frac{x\ln x}{e^x}\,dx=H_1-\frac{G_1}{e^x}+\int\frac x{e^x}\,dx$$ Similarly, if we denote $I_n$ as the integral of $xe^{-x}\ln^nx$ and $H_n$ as the integral of $e^{-x}\ln^nx$, we get $$I_2=H_2-\frac{G_2}{e^x}+2I_1-2\int\frac x{e^x}\,dx$$ and $$I_3=H_3-\frac{G_3}{e^x}+3I_2-6I_1+6\int\frac x{e^x}\,dx$$ which can be written as $$I_3=H_3+3H_2-\frac1{e^x}(G_3-3G_2)$$
Step III:
Now integrate by parts $I_3$. We will integrate the polynomial part of the integrand ($x$) and differentiate the rest ($e^{-x}\ln^3x$). So $$I_3=\left[\frac{x^2\ln^3x}{2e^x}\right]-\int\frac{x\ln^2x(3-x\ln x)}{2e^x}\,dx=\frac{x^2\ln^3x}{2e^x}-\frac32I_2+\frac12\int\frac{x^2\ln^3x}{e^x}\,dx$$ giving $$J_3=\int\frac{x^2\ln^3x}{e^x}\,dx=2I_3+3I_2-\frac{x^2\ln^3x}{e^x}$$ Hence your indefinite integral is $$K=\int\frac{(x^2-3x+1)\ln^3x}{e^x}\,dx=J_3-3I_3+H_3=\frac1{e^x}(G_3-3G_2-x^2\ln^3x)+3H_2$$
A: 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. $$
