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.
Probably not the most direct way to obtain the result, but one possible method:
From the integral representation DLMF
\begin{equation}
\int_{0}^{\infty}\frac{e^{-a^2t}}{\sqrt{t+\cos^{2}\theta}}\mathrm{d}t=\frac{\sqrt{\pi}}{a}e^{a^2\cos^{2}\theta}\operatorname{erfc}\left(a\cos\theta\right) \tag{1}\label{1}
\end{equation}
the integral can be transformed into
\begin{equation}
I=\frac{a}{\sqrt{\pi}}\int_{0}^{\pi/2}\cos \theta \,d\theta \int_{0}^{\infty}\frac{e^{-a^2t}}{\sqrt{t+\cos^{2}\theta}}\,\mathrm{d}t
\end{equation}
By changing the integration order,
\begin{align}
I&=\frac{a}{\sqrt{\pi}}\int_{0}^{\infty}e^{-a^2t}\,dt\int_{0}^{\pi/2}\frac{\cos \theta }{\sqrt{t+1-\sin^2\theta}}\,d\theta\\
&=\frac{a}{\sqrt{\pi}}\int_{0}^{\infty}e^{-a^2t}\arcsin \frac{1}{\sqrt{t+1}}\,dt
\end{align}
Now, integrating by parts, ($u'=e^{-a^2t},v=\arcsin \frac{1}{\sqrt{t+1}}$), we get
\begin{align}
I&=\frac{a}{\sqrt{\pi}}\left[\frac{\pi}{2a^2}-\frac{1}{2a^2}\int_{0}^{\infty}\frac{e^{-a^2t}}{\sqrt{t}(t+1)} \right]\\
&=\frac{a}{\sqrt{\pi}}\left[\frac{\pi}{2a^2}-\frac{1}{a^2}\int_{0}^{\infty}\frac{e^{-a^2u^2}}{u^2+1} \right]
\end{align}
(by changing $t=u^2$). Using the integral representation DLMF
\begin{equation}
\operatorname{erfc}(a)=\frac{2}{\pi}e^{-a^{2}}\int_{0}^{\infty}\frac{e^{-a^{2}u^%
{2}}}{u^{2}+1}\mathrm{d}u \tag{2}\label{2}
\end{equation}
the final result is obtained:
\begin{equation}
I= \frac{\sqrt{\pi } \left(1-e^{a^2} \operatorname{erfc}\left(a\right)\right)}{2 a}
\end{equation}
\eqref{1} This identity is retrieved by changing $t=u^2-\cos^2\theta$ in the integral.
\eqref{2}
With $A=a^2$, defining
\begin{align}
J(A)&=\frac{2}{\pi}\int_0^\infty\frac{e^{-A\left( u^2+1 \right)}}{u^2+1}\\
\frac{dJ(A)}{dA}&=-\frac{2}{\pi}\int_0^\infty e^{-A\left( u^2+1 \right)}\,du\\
&=-\frac{1}{\sqrt{\pi A}}e^{-A}
\end{align}
and remarking that $J(0)=1$,
\begin{align}
J(A)&=1-\frac{1}{\sqrt{\pi}}\int_0^A\frac{e^{-s}}{\sqrt{s}}\,ds\\
&=1-\operatorname{erf}(a)
\end{align}
Best Answer
If you substitute $av=x$ you are left with the problem of evaluating an integral of the form $$ \int_z^{ + \infty } {\frac{\operatorname{erfc}x}{{x^2 }}dx},\quad z>0 . $$ Using integration by parts once yields $$ \int_z^{ + \infty } {\frac{{\operatorname{erfc}x}}{{x^2 }}dx} = \frac{{\operatorname{erfc}z}}{z} - \frac{2}{{\sqrt \pi }}\int_z^{ + \infty } {\frac{1}{x}e^{ - x^2 } dx} = \frac{{\operatorname{erfc}z}}{z} - \frac{1}{{\sqrt \pi }}E_1 (z^2 ), $$ where $E_1$ is the exponential integral.