Let $u=\text{erf}^{-1}(x)$ so that $x=\text{erf}\left(u\right)$, and $dx=\frac{2}{\sqrt{\pi}}e^{-u^{2}}du.$ Then $$I(\sigma)=\frac{1}{2}\int_{-1}^{1}\text{erf}\left(\frac{\sigma}{\sqrt{2}}+\text{erf}^{-1}(x)\right)dx=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-u^{2}}\text{erf}\left(\frac{\sigma}{\sqrt{2}}+u\right)du.$$ Differentiating with respect to $\sigma$ we find that $$\frac{\partial I(\sigma)}{\partial\sigma}=\frac{\sqrt{2}}{\pi}\int_{-\infty}^{\infty}e^{-u^{2}}e^{-(u+\sigma/\sqrt{2})^{2}}du,$$ and so substituting $u=y/\sqrt{2}$ we have $$\frac{\partial I(\sigma)}{\partial\sigma}=\frac{1}{\pi}\int_{-\infty}^{\infty}e^{-y^{2}+y\sigma-\sigma^{2}/2}dy=\frac{1}{\pi}\int_{-\infty}^{\infty}e^{-(y+\sigma/2)^{2}-\sigma^{2}/4}dy=\frac{1}{\sqrt{\pi}}e^{-\sigma^{2}/4}.$$ Hence $$I(\sigma)=\frac{1}{\sqrt{\pi}}\int_{0}^{\sigma}e^{-t^{2}/4}dt=\frac{2}{\sqrt{\pi}}\int_0^{\sigma/2}e^{-t^2}dt=\text{erf}\left(\sigma/2\right),$$ as desired.
The result reads:
\begin{eqnarray}
&&\int\limits_0^\infty x e^{-a x^2-b x} \text{erf}(c x+d) dx= \frac{e^{\frac{b^2}{4 a}}}{4 a} \left( \right.\\
&& \frac{2 b \left(\arctan\left(\frac{\sqrt{a} (b+2 c d)}{2 a d-b c}\right)+\arctan\left(\frac{2 \sqrt{a} d}{b}\right)-\arctan\left(\frac{c}{\sqrt{a}}\right)\right)}{\sqrt{\pi } \sqrt{a}} +\\
&& 2 e^{-\frac{b^2}{4 a}} \text{erf}(d)-\frac{\sqrt{\pi } b \text{erf}\left(\frac{2 a d-b c}{2 \sqrt{a \left(a+c^2\right)}}\right)}{\sqrt{a}}+\frac{2 c e^{-\frac{(b c-2 a d)^2}{4 a \left(a+c^2\right)}}
\text{erfc}\left(\frac{b+2 c d}{2 \sqrt{a+c^2}}\right)}{\sqrt{a+c^2}} +\\
&& -\frac{4 \sqrt{\pi } b \left(T\left(\frac{2 a d-b c}{\sqrt{2} \sqrt{a \left(a+c^2\right)}},\frac{\sqrt{a} (b+2 c d)}{2 a d-b c}\right)+T\left(\frac{b}{\sqrt{2} \sqrt{a}},\frac{2 \sqrt{a}
d}{b}\right)\right)}{\sqrt{a}}\\
&&\left.\right)
\end{eqnarray}
where $T(\cdot,\cdot)$ is the Owen's t function https://en.wikipedia.org/wiki/Owen%27s_T_function .
For[count = 1, count <= 200, count++,
{a, b, c} = RandomReal[{0, 5}, 3, WorkingPrecision -> 50];
d = RandomReal[{-5, 5}, WorkingPrecision -> 50];
I1 = NIntegrate[x Exp[-a x^2 - b x] Erf[c x + d], {x, 0, Infinity},
WorkingPrecision -> 20];
1/a NIntegrate[
x Exp[-x^2 - b/Sqrt[a] x] Erf[c/Sqrt[a] x + d], {x, 0, Infinity}];
1/a Exp[
b^2/(4 a)] NIntegrate[(x - b/(2 Sqrt[a])) Exp[-x^2 ] Erf[
c/Sqrt[a] ( x - b/(2 Sqrt[a])) + d], {x, b/(2 Sqrt[a]),
Infinity}];
1/a Exp[
b^2/(4 a)] (NIntegrate[(x) Exp[-x^2 ] Erf[
c/Sqrt[a] ( x - b/(2 Sqrt[a])) + d], {x, b/(2 Sqrt[a]),
Infinity}] -
b Sqrt[Pi]/(2 Sqrt[a]) 2 T[b/Sqrt[2 a],
c/Sqrt[a], (-((b c)/(2 a)) + d) Sqrt[2]]);
I2 = E^(b^2/(4 a))/(
4 a) ((2 b)/(
Sqrt[a] Sqrt[\[Pi]]) (-ArcTan[c/Sqrt[a]] +
ArcTan[(2 Sqrt[a] d)/b] +
ArcTan[(Sqrt[a] (b + 2 c d))/(-b c + 2 a d)]) +
2 E^(-(b^2/(4 a))) Erf[d] - (b Sqrt[\[Pi]])/Sqrt[a]
Erf[(-b c + 2 a d)/(2 Sqrt[a (a + c^2)])] + (2 c)/Sqrt[
a + c^2] E^(-((b c - 2 a d)^2/(4 a (a + c^2))))
Erfc[(b + 2 c d)/(2 Sqrt[a + c^2])] - (4 b Sqrt[\[Pi]])/Sqrt[
a] (OwenT[b/(Sqrt[2] Sqrt[a]), (2 Sqrt[a] d)/b] +
OwenT[(-b c + 2 a d)/(Sqrt[2] Sqrt[a (a + c^2)]), (
Sqrt[a] (b + 2 c d))/(-b c + 2 a d)]));
If[Abs[I2/I1 - 1] > 10^(-3),
Print["results do not match", {a1, a2, c, d, {I1, I2}}]; Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
];
Best Answer
Your integral admits a closed form.
Theorem. Let $a$ and $b$ be any real numbers. Then
where $\displaystyle \text{erf}(z)=\frac{2}{\sqrt{\pi }}\int _0^ze^{\large-t^2}\mathrm{d}t.$
Proof. Let's denote the left hand side of $(1)$ by $I(a,b)$.
By the change of variable $$\displaystyle U:=\frac{a-\log x}{b}, \quad dx=-b\:e^a \:e^{-b U}dU,$$ we have $$ \begin{align} I(a,b)&=|b|\:e^a \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \left(\text{erf}\left(U\right)+1\right) \mathrm{d}U\\\\ &=|b|\:e^a \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \text{erf}(U)\: \mathrm{d}U+|b|\:e^a \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \mathrm{d}U \tag2 \end{align} $$ Clearly, by the gaussian integral: $$ \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \mathrm{d}U =e^{\large b^2/4}\int_{-\infty }^{+\infty } e^{\large-(U+b/2)^2} \mathrm{d}U =\sqrt{\pi}\:e^{ b^2/4} .\tag3 $$ Set $$ f(b):=\int_{-\infty }^{+\infty } e^{\large-U^2-bU} \text{erf}(U)\: \mathrm{d}U \tag4 $$ and observe that $$ f(0)=\frac{\sqrt{\pi}}{4}\left[(\text{erf}(U))^2\right]_{-\infty }^{+\infty }=0. \tag5$$ Differentiating $(4)$ with respect to $b$ and performing an integration by parts gives $$ \begin{align} f'(b)&=-\int_{-\infty }^{+\infty } U\:e^{\large-U^2-bU} \text{erf}(U)\:\mathrm{d}U\\\\ f'(b)&=\left[\frac 12 e^{\large-U^2}\left(e^{\large-bU} \text{erf}(U)\right)\right]_{-\infty }^{+\infty }-\frac 12\int_{-\infty }^{+\infty } e^{\large-U^2} \left(-b\:e^{\large-bU} \text{erf}(U)+\frac{2}{\sqrt{\pi }}e^{\large-bU}e^{\large-U^2}\right)\mathrm{d}U\\\\ f'(b)&=\frac{b}{2}f(b)-\frac{1}{\sqrt{\pi }}\int_{-\infty }^{+\infty } e^{\large-2U^2-bU}\mathrm{d}U \end{align} $$ or $$ f'(b)=\frac{b}{2}f(b)-\frac{\sqrt{2}}{2}\:e^{\large b^2/8}. \tag6 $$ We classically solve the ODE $(6)$, using $(5)$, to obtain $$ f(b)=-\sqrt{\pi }\:e^{\large b^2/4}\text{erf}\!\left(\! \frac{\sqrt{2}\:b}{4}\!\right) \tag7 $$ then plugging $(7)$, $(4)$ and $(3)$ into $(2)$ gives $(1)$ as desired.