Closed form for $\int_{0}^{\infty }\!{\rm erf} \left(cx\right) \left( {\rm erf} \left(x \right) \right) ^{2}{{\rm e}^{-{x}^{2}}}\,{\rm d}x$
I think it is more practical to consider the parametric integral $$ I(a,b,c) = \int_{0}^{+\infty}\!\!\!\text{erf}(ax)\,\text{erf}(bx)\,\text{erf}(cx)\,e^{-x^2}dx \tag{1}$$ and notice that: $$ \frac{\partial^3 I}{\partial a\,\partial b\,\partial c}=\frac{8}{\pi^{3/2}}\int_{0}^{+\infty}x^3 e^{-(1+a^2+b^2+c^2)x^2}\,dx=\frac{4}{\pi^{3/2}(1+a^2+b^2+c^2)^2}\tag{2}$$ Integrating over $(b,c)\in(0,1)^2$ then with respect to $a$ is now tedious but doable with the support of a CAS. Obviously, we would be happier to integrate over $0\leq b^2+c^2\leq 2$ or $0\leq b^2+c^2\leq 1$ for first: this change of integration domain leads to decent bounds for $I(a,b,c)$. Anyway, the problem is equivalent to integrating a Cauchy-like distribution over $R=[0,1]\times[0,1]\times[0,a]$: a probabilistic approach through characteristic functions might be interesting, too.
Partial work. Consider $$I\left(c\right)=\int_{0}^{\infty}\textrm{erf}\left(cx\right)\textrm{erf}^{2}\left(x\right)e^{-x^{2}}dx. $$ Then we have, integrating by parts, that $$I\left(c\right)=\frac{\sqrt{\pi}}{6}-\frac{c}{3}\int_{0}^{\infty}\textrm{erf}^{3}\left(x\right)e^{-c^{2}x^{2}}dx=\frac{\sqrt{\pi}}{6}-\frac{c}{3}\int_{0}^{\infty}\left(1-\textrm{erfc}\left(x\right)\right)^{3}e^{-c^{2}x^{2}}dx $$ $$=\frac{\sqrt{\pi}}{6}-\frac{c}{3}\int_{0}^{\infty}e^{-c^{2}x^{2}}dx+\frac{c}{3}\int_{0}^{\infty}\textrm{erfc}^{3}\left(x\right)e^{-c^{2}x^{2}}dx $$ $$-c\int_{0}^{\infty}\textrm{erfc}^{2}\left(x\right)e^{-c^{2}x^{2}}dx+c\int_{0}^{\infty}\textrm{erfc}\left(x\right)e^{-c^{2}x^{2}}dx $$ $$=\frac{\sqrt{\pi}}{6}-\frac{c}{3}H_{0,0}\left(c^{2}\right)+\frac{c}{3}H_{0,3}\left(c^{2}\right)-cH_{0,2}\left(c^{2}\right)+cH_{0,1}\left(c^{2}\right), $$ say. I used this particular notation because now we may use these results at page $13$ and get $$H_{0,0}\left(c^{2}\right)=\frac{\sqrt{\pi}}{2c} $$ $$H_{0,1}\left(c^{2}\right)=\frac{\arctan\left(c\right)}{c\sqrt{\pi}} $$ $$H_{0,2}\left(c^{2}\right)=\frac{1}{c\sqrt{\pi}}\left(2\arctan\left(c\right)-\arccos\left(\frac{1}{1+c^{2}}\right)\right) $$ and, using the corollary $10.7$, $$H_{0,3}\left(c^{2}\right)=\frac{\sqrt{\pi}}{2c}-\frac{\arctan\left(1/c\right)}{c\sqrt{\pi}}+\frac{6}{c\pi}\int_{c^{2}}^{\infty}\frac{\arctan\left(\sqrt{t+2}\right)}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dt. $$ At this moment I don't know if it is possible to evaluate the last integral, I will work on it later. However it is quite simple to prove that $$\int_{c^{2}}^{\infty}\frac{1}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dx=\frac{\pi}{2}-2\arctan\left(\frac{c}{\sqrt{c^{2}+2}}\right) $$ so obviously $$\int_{c^{2}}^{\infty}\frac{\arctan\left(\sqrt{t+2}\right)}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dt\leq\frac{\pi}{4}-\pi\arctan\left(\frac{c}{\sqrt{c^{2}+2}}\right) $$ and $$\int_{c^{2}}^{\infty}\frac{\arctan\left(\sqrt{t+2}\right)}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dt\geq\left(\frac{\pi}{2}-2\arctan\left(\frac{c}{\sqrt{c^{2}+2}}\right)\right)\arctan\left(\sqrt{c^{2}+2}\right) $$ so it is not a closed form but it could help.
Let us denote : \begin{equation} {\mathcal I}(c) := \int\limits_0^\infty \operatorname{erf}(c x) \cdot [\operatorname{erf}( x)]^2 e^{-x^2} dx \end{equation} By differentiating with respect to the parameter $c$ we have: \begin{equation} \frac{ d }{d c} {\mathcal I}(c) = \frac{2^2}{\pi^{3/2}} \frac{1}{1+c^2} \cdot \frac{1}{\sqrt{2+c^2}} \cdot \arctan\left( \frac{1}{\sqrt{2+c^2}}\right) \end{equation} therefore the only thing we need to do is to integrate the right hand side. I have calculated a more generic integral that involves this one as a special case in A generalized Ahmed's integral . Here I only state the result: \begin{eqnarray} &&{\mathcal I}(c) = \frac{4}{\pi^{3/2}} \left(\right.\\ && \arctan( \frac{c}{\sqrt{2+c^2}}) \arctan( \frac{1}{\sqrt{2+c^2}})+\\ && \frac{\imath}{2} \left.\left[ {\mathcal F}^{(\alpha_-,+e^{-\imath \phi})}(t)+ {\mathcal F}^{(\alpha_-,-e^{-\imath \phi})}(t)- {\mathcal F}^{(\alpha_-,-e^{+\imath \phi})}(t)- {\mathcal F}^{(\alpha_-,+e^{+\imath \phi})}(t) \right]\right|_0^B-\\ && \frac{\imath}{2} \left.\left[ {\mathcal F}^{(\alpha_+,+e^{-\imath \phi})}(t)+ {\mathcal F}^{(\alpha_+,-e^{-\imath \phi})}(t)- {\mathcal F}^{(\alpha_+,-e^{+\imath \phi})}(t)- {\mathcal F}^{(\alpha_+,+e^{+\imath \phi})}(t) \right]\right|_0^B \left.\right) \end{eqnarray} where $\alpha_- = \sqrt{2}-1$, $\alpha_+:=\sqrt{2}+1$, $\phi:= \arccos(1/\sqrt{3})$, $B:=(-\sqrt{2}+\sqrt{2+c^2})/c$ and \begin{eqnarray} &&{\mathcal F}^{(a,b)}(t):=\int \arctan(\frac{t}{a}) \frac{1}{t-b} dt = \log(t-b) \arctan(\frac{t}{a})\\ &&-\frac{1}{2 \imath} \left( \log(t-b) \left[ \log(\frac{t-\imath a}{b-\imath a}) - \log(\frac{t+\imath a}{b+\imath a})\right] + Li_2(\frac{b-t}{b-\imath a}) - Li_2(\frac{b-t}{b+\imath a})\right) \end{eqnarray}
Update:
Note that the anti-derivative ${\mathcal F}^{(a,b)}(t)$ may have a jump. This will happen if and only if either the quantity $(t+\imath a)/(b+\imath a)$ crosses the negative real axis or the quantity $(t-\imath a)/(b-\imath a)$ crosses the negative real axis both for some $t\in(0,B)$. This has an effect that the argument of the logarithm jumps by $2\pi$. In order to take this into account we have to exclude from the integration region a small vicinity of the singularity in question. In other words the correct formula reads:
\begin{eqnarray} &&{\mathcal I}(c) = \frac{4}{\pi^{3/2}} \left(\right.\\ && \arctan( \frac{c}{\sqrt{2+c^2}}) \arctan( \frac{1}{\sqrt{2+c^2}})+\\ && \frac{\imath}{2}\left[ {\bar {\mathcal F}}^{(\alpha_-,+e^{-\imath \phi})}(0,B)+ {\bar {\mathcal F}}^{(\alpha_-,-e^{-\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_-,-e^{+\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_-,+e^{+\imath \phi})}(0,B) \right]-\\ && \frac{\imath}{2} \left[ {\bar {\mathcal F}}^{(\alpha_+,+e^{-\imath \phi})}(0,B)+ {\bar {\mathcal F}}^{(\alpha_+,-e^{-\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_+,-e^{+\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_+,+e^{+\imath \phi})}(0,B) \right] \left.\right) \end{eqnarray}
where
\begin{eqnarray} {\bar {\mathcal F}}^{a,b}(0,B) &:=& {\mathcal F}^{(a,b)}(B)-{\mathcal F}^{(a,b)}(A) +\\ && 1_{t^{(*)}_+ \in (0,1)} \left( -{\mathcal F}^{(a,b)}(B(t^{(*)}_+ +\epsilon))+{\mathcal F}^{(a,b)}(B(t^{(*)}_+ -\epsilon))\right)+\\ && 1_{t^{(*)}_- \in (0,1)} \left( -{\mathcal F}^{(a,b)}(B(t^{(*)}_- +\epsilon))+{\mathcal F}^{(a,b)}(B(t^{(*)}_- -\epsilon))\right) \end{eqnarray}
where
\begin{eqnarray} t^{(*)}_\pm:= \frac{Im[\mp \imath a(\bar{b} \mp \imath a)]}{B Im[\bar{b} \mp \imath a]} \end{eqnarray}
See Mathematica code below for testing:
F[t_, a_, b_] :=
Log[t - b] ArcTan[t/a] -
1/(2 I) (Log[
t - b] (Log[(t - I a)/(b - I a)] -
Log[(t + I a)/(b + I a)]) + PolyLog[2, (b - t)/(b - I a)] -
PolyLog[2, (b - t)/(b + I a)]);
FF[A_, B_, a_, b_] := Module[{res, rsp, rsm, tsp, tsm, eps = 10^(-9)},
res = F[B, a, b] - F[A, a, b];
tsp = -(Im[I a (Conjugate[b] - I a)]/(B Im[Conjugate[b] - I a]));
tsm = +(Im[I a (Conjugate[b] + I a)]/(B Im[Conjugate[b] + I a]));
(*
If[0\[LessEqual] tsp\[LessEqual]1,Print["Jump +!!"]];
If[0\[LessEqual] tsm\[LessEqual]1,Print["Jump -!!"]];
*)
rsp = If[
0 <= tsp <= 1, -F[A + (tsp + eps) (B - A), a, b] +
F[A + (tsp - eps) (B - A), a, b], 0];
rsm = If[
0 <= tsm <= 1, -F[A + (tsm + eps) (B - A), a, b] +
F[A + (tsm - eps) (B - A), a, b], 0];
res + rsp + rsm
];
For[count = 1, count <= 100, count++,
c = RandomReal[{-10, 10}, WorkingPrecision -> 50];
x1 = NIntegrate[Erf[c x] Erf[x]^2 Exp[-x^2], {x, 0, Infinity},
WorkingPrecision -> 30];
4/Pi^(3/2)
NIntegrate[
1/(1 + xi^2) 1/Sqrt[2 + xi^2] ArcTan[1/Sqrt[2 + xi^2]], {xi, 0,
c}];
A1 = 1; A2 = 1; A3 = c;
phi = ArcCos[1/Sqrt[3]]; B = (-Sqrt[2] + Sqrt[2 + c^2])/c;
4/Pi^(3/
2) ((ArcTan[c/Sqrt[2 + c^2]]) ArcTan[1 /Sqrt[2 + c^2]] +
4 Sqrt[2]
NIntegrate[(ArcTan[t/(Sqrt[2] - 1)] -
ArcTan[t/(Sqrt[2] + 1)])
t/((1 - t^2)^2 + (2 ) (1 + t^2)^2), {t,
0, (-Sqrt[2] + Sqrt[2 + c^2])/c}, WorkingPrecision -> 30]);
4/Pi^(3/
2) ((ArcTan[c/Sqrt[2 + c^2]]) ArcTan[1 /Sqrt[2 + c^2]] +
I/2 NIntegrate[(ArcTan[t/(Sqrt[2] - 1)] -
ArcTan[t/(Sqrt[2] + 1)]) (1/(t - E^(-I phi)) -
1/(t - E^(I phi)) + 1/(t + E^(-I phi)) -
1/(t + E^(I phi))), {t, 0, (-Sqrt[2] + Sqrt[2 + c^2])/c},
WorkingPrecision -> 30]);
x2 = 4/Pi^(3/2) ((ArcTan[c/Sqrt[2 + c^2]]) ArcTan[1/Sqrt[2 + c^2]] +
I/2 (FF[0, B, (Sqrt[2] - 1), 1/Sqrt[3] - I Sqrt[2/3]] +
FF[0, B, (Sqrt[2] - 1), -(1/Sqrt[3]) + I Sqrt[2/3]] -
FF[0, B, (Sqrt[2] - 1), -(1/Sqrt[3]) - I Sqrt[2/3]] -
FF[0, B, (Sqrt[2] - 1), 1/Sqrt[3] + I Sqrt[2/3]]) -
I/2 (FF[0, B, (Sqrt[2] + 1), 1/Sqrt[3] - I Sqrt[2/3]] +
FF[0, B, (Sqrt[2] + 1), -(1/Sqrt[3]) + I Sqrt[2/3]] -
FF[0, B, (Sqrt[2] + 1), -(1/Sqrt[3]) - I Sqrt[2/3]] -
FF[0, B, (Sqrt[2] + 1), 1/Sqrt[3] + I Sqrt[2/3]]));
If[Abs[x2/x1 - 1] > 10^(-3),
Print["results do not match..", {c, {x1, x2}}]; Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
];