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]];
];
