With some help from Mathematica I got this result:
$$\int_0^\infty\left(5\,x^5+x\right)\operatorname{erfc}\left(x^5+x\right)dx=J_{\frac25}\left(\frac8{25\sqrt5}\right)\left(\frac8{375}\sqrt{\frac25\left(5-\sqrt5\right)}\,\pi\,J_{\frac45}\left(\frac8{25\sqrt5}\right)-\frac8{375}\sqrt{\frac25\left(5+\sqrt5\right)}\,\pi\,J_{\frac65}\left(\frac8{25\sqrt5}\right)\right)+J_{-\frac25}\left(\frac8{25\sqrt5}\right)\left(\frac1{75}\sqrt{2\left(5+\sqrt5\right)}\,\pi\,J_{-\frac15}\left(\frac8{25\sqrt5}\right)+\frac1{25}\left(\sqrt5-5\right)\sqrt{\frac1{10}\left(5+\sqrt5\right)}\,\pi\,J_{\frac15}\left(\frac8{25\sqrt5}\right)-\frac8{375}\sqrt{\frac25\left(5+\sqrt5\right)}\,\pi\,J_{\frac45}\left(\frac8{25\sqrt5}\right)-\frac{4\left(\sqrt5-5\right)\sqrt{2\left(5+\sqrt5\right)}\,\pi\,J_{\frac65}\left(\frac8{25\sqrt5}\right)}{1875}\right)+J_{-\frac15}\left(\frac8{25\sqrt5}\right)\left(\frac1{25}\sqrt{2\left(5-\sqrt5\right)}\,\pi\,J_{\frac25}\left(\frac8{25\sqrt5}\right)+\frac8{375}\sqrt{\frac25\left(5+\sqrt5\right)}\,\pi\,J_{\frac35}\left(\frac8{25\sqrt5}\right)-\frac8{375}\sqrt{\frac25\left(5-\sqrt5\right)}\,\pi\,J_{\frac75}\left(\frac8{25\sqrt5}\right)\right)+J_{\frac15}\left(\frac8{25\sqrt5}\right)\left(-\frac1{75}\sqrt{2\left(5+\sqrt5\right)}\,\pi\,J_{\frac25}\left(\frac8{25\sqrt5}\right)+\frac{4\left(\sqrt5-5\right)\sqrt{2\left(5+\sqrt5\right)}\,\pi\,J_{\frac35}\left(\frac8{25\sqrt5}\right)}{1875}+\frac8{375}\sqrt{\frac25\left(5+\sqrt5\right)}\,\pi\,J_{\frac75}\left(\frac8{25\sqrt5}\right)\right),$$
where $J_\nu(x)$ is the Bessel function of the first kind.
A solution outline:
Note that $(5\,x^5+x)=x\frac{d}{dx}(x^5+x)$. Change the integration variable $y=x^5+x$, then the integral takes the form
$$\int_0^\infty\mathcal{BR}(y)\cdot\operatorname{erfc}y\,dy,$$
where $y\mapsto\mathcal{BR}(y)$ is the inverse (properly selected to satisfy $\mathcal{BR}(y)>0$ for $y>0$) of the polynomial function $x \mapsto x^5+x$. This is a well-known non-elementary function called Bring radical, it can be used to express solutions of quintic equations in an explicit form.
An important fact, it has a representation via a generalized hypergeometric function (I used it to answer another question awhile ago). If we plug the hypergeometric representation into the integral and feed it to Mathematica, it produces the result in terms of Bessel functions shown above. I leave this step without a rigorous proof and rely on Mathematica here. The result agrees with a numerical integration to a very high precision. I would be very glad if anybody could write down an explicit derivation of the formula.
Here is my trial, which is partially successful but still not fully answering to your question.
Using some coordinate change, I derived that
$$ I_{n} := \int_{0}^{\infty} \mathrm{erfc}^{n}(x) \, dx = \frac{1}{\sqrt{n}} \left( \frac{2}{\sqrt{\pi}} \right)^{n-1} \int_{T^{n-1}} \int_{0}^{\infty} s^{n-1}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds d\sigma_{x}, \tag{*} $$
where $T^{n-1} = \{ x \in [0, \infty)^{n} : x_{1} + \cdots + x_{n} = \sqrt{n} \}$ is an $(n-1)$-simplex and $d\sigma_{x}$ is the surface measure on $T^{n-1}$.
Case $n = 1$. Since $T^{0} = \{1\}$ and $d\sigma_{x} = \delta(x-1)$ is a unit mass located at $x = 1$, the equation $\text{(*)}$ gives no new information.
Case $n = 2.$ It is easy to find that
$$ \int_{0}^{\infty} s e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2|x|(|x|+1)}$$
and by suitable parametrization of the line segment $T^{1}$, we get
\begin{align*}
I_{2}
&= \sqrt{\frac{2}{\pi}} \int_{T^{1}} \frac{1}{2|x|(|x|+1)} \, d\sigma_{x} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2}} \frac{\sqrt{2} \, dt}{2\sqrt{t^{2} + (\sqrt{2}-t)^{2}} ( \sqrt{t^{2} + (\sqrt{2}-t)^{2}} + 1)} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\frac{\pi}{4}} \frac{d\theta}{1+\cos\theta}
= \sqrt{\frac{2}{\pi}} \tan\left(\frac{\pi}{8}\right)
= \frac{2 - \sqrt{2}}{\sqrt{\pi}}.
\end{align*}
We can also write
$$ I_{2} = \boxed{ \displaystyle \frac{2}{\sqrt{\pi}} - \frac{2\sqrt{2}}{\pi^{3/2}} \arctan(\infty) }. $$
Case $n = 3$. We have
$$ \int_{0}^{\infty} s^{2}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2\sqrt{\pi}(|x|^{2}-1)} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right). $$
Now we know that $T^{2}$ is a regular triangle with side length $\sqrt{6}$. Thus by introducing polar coordinate change, we get
\begin{align*}
I_{3}
&= \frac{2}{\pi^{3/2}\sqrt{3}} \int_{T^{2}} \frac{1}{|x|^{2}-1} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right) \, d\sigma_{x} \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \int_{0}^{\frac{\sec\theta}{\sqrt{2}}} \frac{1}{r} \left( \frac{\arctan r}{r} - \frac{1}{r^{2} + 1} \right) \, drd\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \left( 1- \frac{\arctan\left(\sec\theta / \sqrt{2}\right)}{\sec\theta / \sqrt{2}} \right) \, d\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}}.
\end{align*}
I haven't tried the last integral, but Mathematica suggests that
$$ \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}} = \frac{\sqrt{3}\pi}{4} -\sqrt{\frac{3}{2}} \arctan\left( \sqrt{2} \right).$$
This gives
$$ I_{3} = \boxed{ \displaystyle \frac{3}{\sqrt{\pi}} - \frac{6\sqrt{2}}{\pi^{3/2}} \arctan\left( \sqrt{2} \right) }, $$
which can be numerically checked.
I am trying to generalize this calculation for higher dimensions, but it seems somewhat daunting. For example, if $n = 4$ we have to evaluate
$$ \int_{0}^{\infty} \mathrm{erfc}^{4}(x) \, dx = \frac{1}{\pi^{3/2}} \int_{T^{3}} \frac{2|x|+1}{(|x|+1)^{2}|x|^{3}} \, d\sigma_{x}. $$
Nevertheless, using the formula
$$ \int_{T^{n-1}} f(|x|) \, d\sigma_{x} = n! \int_{0}^{\sqrt{\frac{1}{n-1}}} \int_{0}^{\sqrt{\frac{n}{n-2}}s_{1}} \cdots \int_{0}^{\sqrt{\frac{3}{1}}s_{n-2}} f\left( \sqrt{1+s_{1}^{2} + \cdots + s_{n-1}^{2}} \right) \, ds_{n-1} \cdots ds_{1} $$
together with aid of Mathematica, I was able to evaluate $I_{4}$ and it was
$$ I_{4} = \boxed{ \displaystyle \frac{4}{\sqrt{\pi}} - \frac{24\sqrt{2}}{\pi^{3/2}} \arctan \left( \frac{1}{\sqrt{8}} \right) }. $$
These series of observations lead us to believe that $I_{n}$ is of the form
$$ I_{n} = \frac{n}{\sqrt{\pi}} - \frac{n!\sqrt{2}}{\pi^{3/2}} \arctan \alpha_{n} $$
for some reasonably nice $\alpha_{n}$ (with $\alpha_{1} = 0$, $\alpha_{2} = \infty$, $\alpha_{3} = \sqrt{2}$ and $\alpha_{4} = \frac{1}{\sqrt{8}}$), but inverse symbolic calculations show that this seems no longer the case for $n \geq 5$.
Indeed, for $n = 5$ we have
$$ I_{5} = \frac{5}{\sqrt{\pi}} - \frac{80\sqrt{2}}{\pi^{3/2}} \arctan \sqrt{\frac{2}{3}} + \frac{240\sqrt{2}}{\pi^{5/2}} A\left( \sqrt{\frac{5}{3}}, \sqrt{\frac{3}{2}}, \sqrt{\frac{4}{5}} \right), $$
where $A(p, q, r)$ is the generalized Ahmed's integral. I have no idea whether it will reduce to a familiar closed form or not.
Best Answer
With the following identity (the proof is easy) $$\mathrm{e}^{x^2}\mathrm{erfc}(x) = \frac{1}{\sqrt{\pi}}\int_0^\infty \mathrm{e}^{-t^2/4}\mathrm{e}^{-xt}\mathrm{d} t, \quad x \ge 0,$$ we have \begin{align} \int_0^\infty \mathrm{e}^{x^2-x} \mathrm{erfc}(x) \mathrm{d} x &= \int_0^\infty \mathrm{e}^{-x} \frac{1}{\sqrt{\pi}}\int_0^\infty \mathrm{e}^{-t^2/4}\mathrm{e}^{-xt}\mathrm{d} t \mathrm{d} x\\ &= \frac{1}{\sqrt{\pi}}\int_0^\infty \mathrm{e}^{-t^2/4}\int_0^\infty \mathrm{e}^{-(1+t)x}\mathrm{d}x \, \mathrm{d} t\\ &= \frac{1}{\sqrt{\pi}}\int_0^\infty \mathrm{e}^{-t^2/4}\frac{1}{1+t} \mathrm{d} t. \end{align} It can be expressed in terms of special functions: $$I = \frac{\pi \mathrm{erfi}(\tfrac{1}{2}) - \mathrm{Ei}(\tfrac{1}{4})}{2\sqrt{\pi}\sqrt[4]{\mathrm{e}}}.$$ Also, it can be expressed in terms of infinite series (for real $x$): $$\mathrm{erfi}(x) = \frac{2}{\sqrt{\pi}}\int_0^x \mathrm{e}^{t^2}\mathrm{d} t = \frac{2}{\sqrt{\pi}}\sum_{k=0}^\infty \frac{x^{2k+1}}{k!(2k+1)}$$ and $$\mathrm{Ei}(x) = \sum_{k=1}^\infty \frac{x^k}{k k!} + \ln x + \gamma.$$