Let us denote:
\begin{equation}
{\mathfrak I}^{(2)}:=\int\limits_{\mathbb R} x \exp\left( -b^2(x-c)^2\right) \cdot [\text{erf}\left( a(x-d)\right)]^2 dx
\end{equation}
Then we have:
\begin{eqnarray}
{\mathfrak I}^{(2)}&=& \int\limits_{\mathbb R} \phi(x) \text{erf}\left(\frac{a b^2 (c-d)}{a^2+b^2}+\frac{a x}{\sqrt{2 a^2+2 b^2}}\right) dx + \frac{c \sqrt{\pi}}{b}\int\limits_{\mathbb R} \phi(x) [\text{erf}\left(\frac{a x}{\sqrt{2} b}+a (c-d)\right)]^2 dx =\\
&& \text{sign}(b) \left(\right.\\
%
&&\frac{2 a \sqrt{\frac{b^2}{a^2+b^2}} e^{-\frac{a^2 b^2 (c-d)^2}{a^2+b^2}} \text{erf}\left(\frac{a b^2 (c-d)}{\sqrt{a^2+b^2} \sqrt{2 a^2+b^2}}\right)}{b^3}+\\
%
&&-\frac{4 \sqrt{\pi } c T\left(\frac{\sqrt{2} a (c-d)}{\sqrt{\frac{a^2}{b^2}+1}},\frac{1}{\sqrt{\frac{2 a^2}{b^2}+1}}\right)}{b}-\frac{4 \sqrt{\pi } c T\left(\frac{a \sqrt{\frac{2 a^2}{b^2}+2}
(c-d)}{\sqrt{\frac{\left(a^2+b^2\right)^2}{b^4}}},\frac{1}{\sqrt{\frac{2 a^2}{b^2}+1}}\right)}{b}+\frac{\sqrt{\pi } c \text{erf}\left(\frac{a (c-d)}{\sqrt{\frac{a^2}{b^2}+1}}\right)}{b}+\frac{4 c \arctan\left(\frac{1}{\sqrt{\frac{2 a^2}{b^2}+1}}\right)}{\sqrt{\pi } b}+\frac{2 c \arctan\left(\frac{a^2}{b^2 \sqrt{\frac{2 a^2}{b^2}+1}}\right)}{\sqrt{\pi } b}-\frac{\sqrt{\pi } c \text{erf}\left(\frac{a
\sqrt{\frac{a^2}{b^2}+1} (c-d)}{\sqrt{\frac{\left(a^2+b^2\right)^2}{b^4}}}\right)}{b}\\
&&\left.\right)
\end{eqnarray}
In the first line we substituted for $x-c$ and then split the integral into two integrals one that involves $x$ and the other that does not. Subsequently we integrated by parts the one that involves $x$. Now in the second line we just used A definite integral involving a Gaussian and shifted error functions. .
Clear[phi]; phi[x_] := Exp[-1/2 x^2]/Sqrt[2 Pi]; eps = 10^(-9);
For[count = 1, count <= 500, count++,
{a, b, c, d} = RandomReal[{-3, 3}, 4, WorkingPrecision -> 50];
x1 = NIntegrate[
x Exp[-b^2 (x - c)^2] Erf[a (x - d)]^2, {x, -Infinity, Infinity}];
Sqrt[2 Pi]/(2 b^2) NIntegrate[(x) phi[
x] Erf[a (c - d) + (a x)/(Sqrt[2] b)]^2, {x, -Infinity,
Infinity}] +
c Sqrt[2 Pi]/(Sqrt[2] b) NIntegrate[
phi[x] Erf[a (c - d) + (a x)/(Sqrt[2] b)]^2, {x, -Infinity,
Infinity}];
(2 a Sqrt[b^2/(a^2 + b^2)] E^(-((a^2 b^2 (c - d)^2)/(a^2 + b^2))))/
b^3 NIntegrate[
phi[x] (
Erf[ ((a (b^2) )/(a^2 + b^2)) (c - d) +
a/Sqrt[2 a^2 + 2 b^2] x]), {x, -Infinity, Infinity}] + (
c Sqrt[\[Pi]])/
b NIntegrate[
phi[x] Erf[a (c - d) + (a x)/(Sqrt[2] b)]^2, {x, -Infinity,
Infinity}];
(2 a Sqrt[b^2/(a^2 + b^2)] E^(-((a^2 b^2 (c - d)^2)/(a^2 + b^2))))/
b^3 2 (T[eps, a/Sqrt[
a^2 + b^2], ((Sqrt[2] a (b^2) )/(a^2 + b^2)) (c - d)] +
T[eps, -(a/Sqrt[
a^2 + b^2]), ((Sqrt[2] a (b^2) )/(a^2 + b^2)) (c - d)]) + (
c Sqrt[\[Pi]])/
b NIntegrate[
phi[x] Erf[a (c - d) + (a x)/(Sqrt[2] b)]^2, {x, -Infinity,
Infinity}];
x2 = Sign[
b] ((2 a Sqrt[b^2/(a^2 + b^2)]
E^(-((a^2 b^2 (c - d)^2)/(a^2 + b^2))))/
b^3 Erf[(a b^2 (c - d))/(
Sqrt[(a^2 + b^2)] Sqrt[2 a^2 + b^2])] + (c Sqrt[\[Pi]])/
b NIntegrate[
phi[x] (Erf[a (c - d) + (a x)/(Sqrt[2] b)]^2), {x, -Infinity,
Infinity}]);
x2 = Sign[
b] ((2 a Sqrt[b^2/(a^2 + b^2)])/
b^3 E^(-((a^2 b^2 (c - d)^2)/(a^2 + b^2)))
Erf[(a b^2 (c - d))/(Sqrt[(a^2 + b^2)] Sqrt[2 a^2 + b^2])] +
c /(b Sqrt[\[Pi]]) (4 ArcTan[1/Sqrt[1 + (2 a^2)/b^2]] +
2 ArcTan[a^2/(Sqrt[1 + (2 a^2)/b^2] b^2)] + \[Pi] Erf[(
a (c - d))/Sqrt[1 + a^2/b^2]] - \[Pi] Erf[(
a Sqrt[1 + a^2/b^2] (c - d))/Sqrt[(a^2 + b^2)^2/b^4]] -
4 \[Pi] OwenT[(Sqrt[2] a (c - d))/Sqrt[1 + a^2/b^2], 1/Sqrt[
1 + (2 a^2)/b^2]] -
4 \[Pi] OwenT[(a Sqrt[2 + (2 a^2)/b^2] (c - d))/
Sqrt[(a^2 + b^2)^2/b^4], 1/Sqrt[1 + (2 a^2)/b^2]]));
If[Abs[x2/x1 - 1] > 10^(-3),
Print["results do not match..", {a, b, {x1, x2}}]; Break[]];
If[Mod[count, 50] == 0, PrintTemporary[count]];
];