How to prove the result of this definite integral?

Enforce the substitution $x\to x^2$. Then, we have

$$\begin{align} \mathcal{I}&=4\int_0^\infty x^2\log(x)e^{-x}\,dx\\\\ &=4\left.\left(\frac{d}{da}\int_0^\infty x^{2+a}e^{-x}\,dx\right)\right|_{a=0}\\\\ &=4\left.\left(\frac{d}{da}\Gamma(a+3)\right)\right|_{a=0}\\\\ &=4\Gamma(3)\psi(3)\\\\ &=4(2!)(3/2-\gamma)\\\\ &=12-8\gamma \end{align}$$

as was to be shown!

Note the we used (i) the integral representation of the Gamma function

$$\Gamma(x)=\int_0^\infty t^{x-1}e^{-t}\,dt$$

(ii) the relationship between the digamma and Gamma functions

$$\psi(x)=\frac{\Gamma'(x)}{\Gamma(x)}$$

and (iii) the recurrence relationship

$$\psi(x+1)=\psi(x)+\frac1x$$

In addition, we used the special values

$$\Gamma(3)=2!$$

and

$$\psi(1)=-\gamma$$

Hint. One may recall that $$ \int_0^\infty u^{s}e^{-u}\:du=\Gamma(s+1), \quad s>0,\tag1 $$ giving, by differentiating under the integral sign, $$ \begin{align} \int_0^\infty u^{2}\ln u \:e^{-u}\:du&=\left.\left(\Gamma(s+1)\right)'\right|_{s=2}\\\\ &=\left.\left(s(s-1)\Gamma(s-1)\right)'\right|_{s=2}\\\\ &=3+2\Gamma'(1)\\\\ &=3-2\gamma,\tag2 \end{align} $$ where we have used $\Gamma'(1)=-\gamma$.

Then, one may rewrite the initial integral as $$ 2\int_0^\infty \sqrt{x}\:(\ln \sqrt{x}) \:e^{-\sqrt{x}}\:dx, $$ then perform the change of variable $x=u^2$, $dx=2udu$, obtaining $$ \begin{align} \int_0^\infty \sqrt{x}\ln x \:e^{-\sqrt{x}}\:dx&=4\int_0^\infty u^{2}\ln u \:e^{-u}\:du.\tag3 \end{align} $$ Considering $(2)$ and $(3)$ gives the announced result.

Start from the well-known integral

$$-\gamma=\int_0^\infty\exp(-x)\log x\,\mathrm dx$$

A round of integration by parts yields

$$\begin{align*} -\gamma&=\int_0^\infty\exp(-x)\log x\,\mathrm dx\\ &=\int_0^\infty x\exp(-x)\log x\,\mathrm dx-\int_0^\infty x\exp(-x)\,\mathrm dx\\ 1-\gamma&=\int_0^\infty x\exp(-x)\log x\,\mathrm dx \end{align*}$$

A second round gives

$$\begin{align*} 1-\gamma&=\int_0^\infty x\exp(-x)\log x\,\mathrm dx\\ &=\int_0^\infty x\exp(-x)(x-1)(\log x-1)\,\mathrm dx\\ &=\int_0^\infty x\exp(-x)\,\mathrm dx-\int_0^\infty x^2\exp(-x)\,\mathrm dx-\int_0^\infty x\exp(-x)\log x\,\mathrm dx+\int_0^\infty x^2\exp(-x)\log x\,\mathrm dx\\ 3-2\gamma&=\int_0^\infty x^2\exp(-x)\log x\,\mathrm dx \end{align*}$$

where the last integral is the one obtained by Dr. MV after an appropriate substitution. Thus, the original integral is equal to $12-8\gamma$.