Integral $\int_0^\infty x^2\,e^{-x^2}\operatorname{erf}(x)\,\log(x)\,dx$
$\newcommand{\+}{^{\dagger}}% \newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}% \newcommand{\dd}{{\rm d}}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\fermi}{\,{\rm f}}% \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}% \newcommand{\half}{{1 \over 2}}% \newcommand{\ic}{{\rm i}}% \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ket}[1]{\left\vert #1\right\rangle}% \newcommand{\ol}[1]{\overline{#1}}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}% \newcommand{\sech}{\,{\rm sech}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
$\ds{I \equiv \int_{0}^{\infty}x^{2}\expo{-x^{2}}{\rm erf}\pars{x}\ln\pars{x} \,\dd x}$. Let's $\ds{{\cal I}\pars{\mu} \equiv \int_{0}^{\infty}x^{\mu}\expo{-x^{2}}{\rm erf}\pars{x}\,\dd x}$ such that $\ds{I = \lim_{\mu \to 2}\totald{{\cal I}\pars{\mu}}{\mu}}$
Since $\ds{{\rm erf}\pars{x} \stackrel{{\rm def.}}{=}{2 \over \root{\pi}}\int_{0}^{x}\expo{-y^{2}}\,\dd y}$: \begin{align} {\cal I}\pars{\mu}&=\int_{0}^{\infty}x^{\mu}\expo{-x^{2}} {2 \over \root{\pi}}\int_{0}^{\infty}\Theta\pars{x - y}\expo{-y^{2}}\,\dd y\,\dd x \\[3mm]&= {2 \over \root{\pi}}\int_{0}^{\pi/2}\dd\theta\,\cos^{\mu}\pars{\theta} \Theta\pars{\cos\pars{\theta} - \sin\pars{\theta}} \overbrace{\int_{0}^{\infty}\dd r\,r^{\mu + 1}\expo{-r^{2}}} ^{\Gamma\pars{1 + \mu/2}/2} \\[3mm]&={1 \over \root{\pi}}\,\Gamma\pars{1 + {\mu \over 2}} \int_{0}^{\pi/4}\dd\theta\,\cos^{\mu}\pars{\theta} \end{align}
\begin{align} I&=\lim_{\mu \to 2}\totald{{\cal I}\pars{\mu}}{\mu}= {1 \over 2\root{\pi}}\,\overbrace{\Psi\pars{2}}^{\ds{1 - \gamma}}\ \overbrace{\int_{0}^{\pi/4}\dd\theta\,\cos^{2}\pars{\theta}}^{\ds{\pars{\pi + 2}/8}} \\[3mm]&+ {1 \over \root{\pi}}\,\overbrace{\Gamma\pars{2}}^{\ds{1}} \overbrace{% \int_{0}^{\pi/4}\dd\theta\,\cos^{2}\pars{\theta}\ln\pars{\cos\pars{\theta}}} ^{\ds{\braces{4G + \pi - 2\bracks{1 + \ln\pars{2}} - \pi\ln\pars{4}}/16}} \end{align} $\Gamma\pars{z}$ and $\Psi\pars{z}$ are the $\it Gamma$ and $\it Digamma$ functions, respectively. $\gamma$ and $G$ are the $\it Euler-Mascheroni$ and Catalan constants, respectively.
$$ \begin{array}{|l|}\hline \mbox{}\\ \quad{\displaystyle\int_{0}^{\infty}x^{2}\expo{-x^{2}} \,\mathrm{erf}\pars{x}\ln\pars{x}\,\dd x}\quad \\[2mm] = \quad{{\displaystyle\quad\pars{\pi + 2}\pars{1 - \gamma} + 4G + \pi - 2\bracks{1 + \ln\pars{2}} - \pi\ln\pars{4}\quad} \over {\displaystyle 16\root{\pi}}}\quad \\ \mbox{}\\ \hline \end{array} \approx 0.0436462 $$
Below is almost a copypaste of my answer to equivalent question - which was asked one month earlier than this one.
We start with the standard integral representation of the hypergeometric function: \begin{align} _2F_1\left(\frac12,a,\frac32,-1\right)=\frac12\int_0^1\frac{dt}{\sqrt{t}\left(1+t\right)^{a}}. \end{align} Differentiating it w.r.t. $a$, one finds \begin{align} S:=\left[\frac{d}{da} {}_2F_1\left(\frac12,a,\frac32,-1\right)\right]_{a=2}= -\frac12\int_0^1\frac{\ln\left(1+t\right)dt}{\sqrt{t}\left(1+t\right)^{2}}. \end{align} After the change of variables $t=s^2$ the last integral can be expressed in terms of dilogarithms: \begin{align} S&=-\int_0^1\frac{\ln\left(1+s^2\right)ds}{\left(1+s^2\right)^2}=\\&= \frac{\pi}{8}\left[1-3\ln 2 +\ln\left(2+\sqrt{2}\right)\right]-\frac{1+\ln 2}{4}+\Im\left(\operatorname{Li}_2\left(-e^{i\pi/4}\right)-\operatorname{Li}_2\left(1-e^{i\pi/4}\right)\right)=\\ &=\frac{\pi\left(1-2\ln 2\right)}{8}-\frac{1+\ln 2}{4}+\frac12 \Im\operatorname{Li}_2\left(i\right)=\\ &=\frac{G}{2}+\frac{\pi\left(1-2\ln 2\right)}{8}-\frac{1+\ln 2}{4}. \end{align} where $G$ denotes the Catalan's constant. Note that $S$ is precisely what we need to compute.
Explanation: At the first step, the only nontrivial integrals (producing dilogarithms) are $\displaystyle \int\frac{\ln(1\pm i s)}{1\mp is}ds$. The other integrals are elementary. At the last step we used the formula $\operatorname{Li}_2(i)=iG-\frac{\pi^2}{48}$ from here.
$$I=\frac{2-\ln2}{16}\sqrt\pi-\frac{\gamma+\ln2}{16\,\sqrt\pi}(\pi+2)+\frac{G}{4\,\sqrt\pi},$$ where $\gamma$ is the Euler-Mascheroni constant: $$\gamma=-\int_0^1\ln(-\ln x)\,dx$$ and $G$ is the Catalan constant:
$$G=-\int_0^1\frac{\ln x}{x^2+1}dx$$