Interesting, isn't it?

Synopsis : Remove the error function by Integration by Parts , then convert the integral into sum of moments of a standard Gaussian distribution.

---

We first integrate by parts with $u = \text{erf}(\sqrt x)$ , $dv =x^2 e^{-2x^2} \, dx$ . So $v = -\frac14 e^{-2x}(2x^2+2x+1)$ .

By definition, $\displaystyle \text{erf}(x) = \frac2{\sqrt \pi} \int_0^x e^{-t^2} \, dt$ . So by Fundamental Theorem of Calculus , $\frac{d}{dx} \text{erf}(x) = \frac{2}{\sqrt \pi} e^{-x^2}$ . And by Differentiation - Chain Rule , $\displaystyle \frac d{dx} \text{erf}(\sqrt x) = \frac{e^{-x}}{\sqrt \pi \sqrt x}$ . Thus,

$\begin{aligned} & \int_0^\infty & x^2 e^{-2x} \text{erf}(\sqrt x) \, dx \\ &=& \int u \, dv \\ &=& uv - \int v \, du \\ &=& \cancelto{0}{\left[ \frac2{\sqrt \pi} \int_0^x e^{-t^2} \, dt \cdot -\frac14 e^{-2x}(2x^2 + 2x + 1) \right]_{x\to0}^{x\to\infty}} + \frac1{4\sqrt \pi} \int_0^\infty e^{-2x}(2x^2 + 2x + 1) \frac{e^{-x}}{ \sqrt x} \, dx \\ &=& \frac1{4\sqrt \pi} \int_0^\infty \frac{e^{-3x}}{\sqrt x} (2x^2 + 2x + 1) \, dx \quad, \quad \color{#69047E}{\text{Let } y = \sqrt{x} \Rightarrow dx = 2y \, dy} \\ &=& \frac1{2\sqrt \pi} \int_0^\infty e^{-3y^2} (2y^4 + 2y^2 + 1) \, dy \quad, \quad \color{#69047E}{\text{Let } z= \sqrt6 y} \\ &=& \frac1{\sqrt6} \frac1{2\sqrt \pi} \int_0^\infty e^{-z^2/2} \left( \frac{z^4}{18} + \frac{z^2}3 + 1 \right) \, dz \\ &=& \frac1{2\sqrt3} \frac1{\sqrt{2\pi}} \int_0^\infty \left( \frac1{18} z^4 e^{-z^2/2} + \frac13 z^2 e^{-z^2/2} + e^{-z^2 /2} \right) \, dz \\ &=& \frac1{4\sqrt3} \frac1{\sqrt{2\pi}} \int_{-\infty}^\infty \left( \frac1{18} z^4 e^{-z^2/2} + \frac13 z^2 e^{-z^2/2} + e^{-z^2 /2} \right) \, dz \\ &=& \frac1{4\sqrt3} \cdot \left ( \frac1{18} E(X^4) + \frac13 E(X^2) + E(X) \right) \quad, \quad \color{#69047E}{\text{where }X\text{ denote the Gaussian distribution } } \\ &=& \frac1{4\sqrt3} \cdot \left( \frac1{18} \cdot 3!! + \frac13 \cdot 1!! + 1 \right) \\ &=& \boxed{ \dfrac{\sqrt 3}8 } \\ \end{aligned}$

---

Note : The $m^\text{th}$ moment of a Gaussian distribution can be found here .

I know of 2 ways, out of which I will give the solution of one(which is easier and better) and an idea of the other.

$\displaystyle \text{erf}(\sqrt{t}) = \frac{2}{\sqrt{\pi}} \int_{0}^{\sqrt{t}}{e^{-m^2}\, \mathrm{d}m}$

Now we consider our integral as -

$\displaystyle I(n) = \int_{0}^{\infty}{e^{-2t} t^2 \text{erf}(n \sqrt{t})\, \mathrm{d}t}$

and thus we need to find $I(1)$

We would use Leibnitz integral rule

$\displaystyle \frac{\mathrm{d}}{\mathrm{d}x} \left (\int_{a(x)}^{b(x)}f(x,t)\,\mathrm{d}t \right) = f(x,b(x))\cdot b'(x) - f(x,a(x))\cdot a'(x) + \int_{a(x)}^{b(x)} \frac{\partial}{\partial x}f(x,t)\; \mathrm{d}t$

and get that $\displaystyle \frac{\mathrm{d}}{\mathrm{d}n}(\text{erf}(n \sqrt{t})) = \frac{2\sqrt{t} e^{-n^2 t}}{\sqrt{\pi}}$

Hence,

$\displaystyle I'(n) = \int_{0}^{\infty}{e^{-2t} t^2 \frac{2\sqrt{t} e^{-n^2 t}}{\sqrt{\pi}}\, \mathrm{d}t}$

This is a simple Laplace transform,

$\displaystyle I'(n) = \frac{2 \Gamma\left(\frac{7}{2}\right)}{\sqrt{\pi} (n^2+2)^{7/2}}$

$\displaystyle \int_{0}^{1}{I'(n) \, \mathrm{d}n} = \int_{0}^{1}{\frac{2 \Gamma\left(\frac{7}{2}\right)}{\sqrt{\pi} (n^2+2)^{7/2}} \, \mathrm{d}n}$

I will not go into the computation of this integral since it is kinda trivial one.

$\displaystyle I(1) - I(0) = \frac{2}{\sqrt{\pi}} \frac{\Gamma\left(\frac{7}{2}\right)}{10 \sqrt{3}}$

$\displaystyle I(1) = \frac{2}{\sqrt{\pi}} \frac{\frac{15\sqrt{\pi}}{8}}{10 \sqrt{3}}$

$\displaystyle = \boxed{\frac{\sqrt{3}}{8}}$

The other method uses the summation definition of error function and it can take quite a lot of energy even with hypergeometric(if not yours, your computer anyway).