$\operatorname{erf}(x)$ is an odd function, therefore,
$$
\begin{align}
\int_{-\infty}^\infty(\operatorname{erf}(a+x)+\operatorname{erf}(a-x))\;\mathrm{d}x
&=\lim_{L\to\infty}\;\int_{-L}^L(\operatorname{erf}(x+a)-\operatorname{erf}(x-a))\;\mathrm{d}x\\
&=\lim_{L\to\infty}\;\int_{-L+a}^{L+a}\operatorname{erf}(x)\;\mathrm{d}x-\lim_{L\to\infty}\;\int_{-L-a}^{L-a}\operatorname{erf}(x)\;\mathrm{d}x\\
&=\lim_{L\to\infty}\;\int_{L-a}^{L+a}\operatorname{erf}(x)\;\mathrm{d}x-\lim_{L\to\infty}\;\int_{-L-a}^{-L+a}\operatorname{erf}(x)\;\mathrm{d}x\\
&=4a\tag{1}
\end{align}
$$
since $\lim\limits_{x\to\infty}\operatorname{erf}(x)=1$ and $\lim\limits_{x\to-\infty}\operatorname{erf}(x)=-1$.
Furthermore,
$$
\begin{align}
\int_{-\infty}^\infty\left(\operatorname{erf}(a+x)\operatorname{erf}(a-x)+1 \right)\;\mathrm{d}x
&=\int_{-\infty}^\infty(\operatorname{erf}(a+x)+1)(\operatorname{erf}(a-x)+1)\;\mathrm{d}x\\
&-\int_{-\infty}^\infty(\operatorname{erf}(a+x)+\operatorname{erf}(a-x))\;\mathrm{d}x\tag{2}
\end{align}
$$
To evaluate
$$
\begin{align}
\int_{-\infty}^\infty(\operatorname{erf}(a+x)+1)(\operatorname{erf}(a-x)+1)\;\mathrm{d}x
&=\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^{a+x}\int_{-\infty}^{a-x}e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t\;\mathrm{d}x
\end{align}
$$
note that $s\le a+x$ and $t\le a-x$; i.e. $s-a\le x\le a-t$ and $s+t\le2a$. Thus,
$$
\begin{align}
\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^{a+x}\int_{-\infty}^{a-x}e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t\;\mathrm{d}x
&=\frac{4}{\pi}\int\int_{s+t\le2a}\int_{s-a}^{a-t}e^{-s^2-t^2}\;\mathrm{d}x\;\mathrm{d}s\;\mathrm{d}t\\
&=\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^\infty(2a-s-t)_+\;e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t
\end{align}
$$
Change variables: $u=(s+t)/\sqrt{2}$ and $v=(s-t)/\sqrt{2}$ so that $s=(u+v)/\sqrt{2}$ and $t=(u-v)/\sqrt{2}$:
$$
\begin{align}
\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^\infty(2a-s-t)_+\;e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t
&=\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^\infty(2a-\sqrt{2}u)_+\;e^{-u^2-v^2}\;\mathrm{d}u\;\mathrm{d}v\\
&=\frac{4}{\sqrt{\pi}}\int_{-\infty}^{\sqrt{2}a}(2a-\sqrt{2}u)\;e^{-u^2}\;\mathrm{d}u\\
&=4a(\operatorname{erf}(\sqrt{2}a)+1)-\frac{4}{\sqrt{\pi}}\int_{-\infty}^{\sqrt{2}a}\sqrt{2}u\;e^{-u^2}\;\mathrm{d}u\\
&=4a(\operatorname{erf}(\sqrt{2}a)+1)-\frac{2\sqrt{2}}{\sqrt{\pi}}\int_{-\infty}^{\sqrt{2}a}\;e^{-u^2}\;\mathrm{d}u^2\\
&=4a(\operatorname{erf}(\sqrt{2}a)+1)+\frac{2\sqrt{2}}{\sqrt{\pi}}\;e^{-2a^2}
\end{align}
$$
Therefore,
$$
\int_{-\infty}^\infty\left(\operatorname{erf}(a+x)+1\right)\left(\operatorname{erf}(a-x)+1\right)\;\mathrm{d}x
=4a\left(\operatorname{erf}(\sqrt{2}a)+1\right)+\frac{2\sqrt{2}}{\sqrt{\pi}}\;e^{-2a^2}\tag{3}
$$
Thus, the convolution of $\operatorname{erf}(x)+1$ with itself is $2x(\operatorname{erf}(x/\sqrt{2})+1)+\frac{2\sqrt{2}}{\sqrt{\pi}}e^{-x^2/2}$.
Subtract $4a$ from $(3)$ using $(1)$ and $(2)$ to get
$$
\int_{-\infty}^\infty\left(\operatorname{erf}(a+x)\operatorname{erf}(a-x)+1 \right)\;\mathrm{d}x
=4a\operatorname{erf}(\sqrt{2}a)+\frac{2\sqrt{2}}{\sqrt{\pi}}\;e^{-2a^2}\tag{4}
$$
My guess is you want either $(3)$ or $(4)$.
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
We will use the fact that: $\mathrm{erf} (x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^{2}} dt$.
Thus, $$\mathrm{erf}^{2} (x) = \frac{4}{\pi} \int_{0}^{x} \int_{0}^{x} e^{-y^2+z^2}dydz.$$
Replace this expression into the integral you have, and you'll end up with:
$$ \int_{0}^{\infty}\mathrm{erf}^2(x)e^{-(x+b)^2}dx = \frac{4}{\pi}\int_{0}^{\infty}\int_{0}^{x}\int_{0}^{x} e^{-((x+b)^2+y^2+z^2)} dzdydx$$
We'll ignore the $\frac{4}{\pi}$ for now, and just attack the integral. This is a candidate for spherical coordinates. However, this region gets VERY hard to turn into spherical coordinates once you involve the $b$ in the equation. Keep in mind we have to use, instead of the traditional $(x,y,z) = ( \rho \cos \theta \sin\phi, \rho \sin\theta \sin\phi, \rho \cos \phi)$, we have accounted for the $(x+b)$ by using $(x + b,y,z) = ( \rho \cos \theta \sin \phi, \rho \sin \theta \sin \phi, \rho \cos \phi)$. This overall does not change the Jacobian value of $\rho ^2 \sin \phi$, as the $x$ term becomes $\rho \cos \theta \sin \phi - b$, which under all partials, the $-b$ goes to $0$. HOWEVER, the region you're bounding now becomes very complicated. The region that this triple integral encloses is simple enough for $\rho$ and $\theta$. However, because the solid is bound by the plane $z = x + y$, we need to make our $\phi$ value $0 \leq \phi \leq z = x + y$. Replacing $z = x + y$ in for spherical coordinates reveals that $\cos \phi = (\sin\theta + \cos\theta) \sin\phi - b.$ To solve for $\phi$, technically, the solution is:
$$\phi = 2(\arctan(\frac{\sqrt{2 - b^2 + 2\cos\theta \sin\theta}}{b - 1}))$$ if $b \neq 1$ and $b \neq (\sin\theta + \cos \theta)\sqrt{2+2\sin\theta \cos\theta} + 2 + 2\sin\theta \cos\theta.$ We are going to assume this is true, because I'm already getting a headache, haha. Nonetheless, this integral will become:
$$\int_{0}^{\infty} \int_{0}^{\frac{\pi}{4}} \int_{0}^{2(\arctan(\frac{\sqrt{2 - b^2 + 2\cos\theta \sin\theta}}{b - 1}))} e^{-\rho ^2} \rho^2 \sin(\phi)d\phi d\theta d\rho.$$
This actually looks worse than it is, because the $\rho$ integral you can take out, as we know its value. $\int_{0}^{\infty} \rho ^2 e^{-\rho ^2} d \rho = \frac{\sqrt{\pi}}{4}.$ Our integral is just in terms of $\phi$ and $\theta.$ Namely:
$$ \frac{\sqrt{\pi}}{4} \int_{0}^{\frac{\pi}{4}} \int_{0}^{2(\arctan(\frac{\sqrt{2 - b^2 + 2\cos\theta \sin\theta}}{b - 1}))} \sin(\phi)d\phi d\theta.$$
This integral actually has an incredibly simple integrand, but it's just the bounds that are confusing. Let's rename $\phi$'s upper bound to $f(\theta)$. Integrating with respect to $\phi$ first, we get:
$$\int_{0}^{\frac{\pi}{4}} [ -\cos(\phi)]_{0,f(\theta)} d\theta = -\int_{0}^{\frac{\pi}{4}} 1 - \cos(f(\theta)) d\theta.$$
Now, this is about where I stop. You can replace $f(\theta)$ with your arctangent business, and take the cosine of the $2 \times$arctangent using some triangles, and I'm sure you will get something nice and neat out of it. I can already see using the double angle formula for cosine being really helpful to cancel square roots. That $\phi$ bound is an absolute headache. Replace what I gave you for $f(\theta)$ and substitute it in, then integrate one more time from $0$ to $\frac{\pi}{4}$ and then multiply by the $\frac{4}{\pi}$ that we got from the start. This integral I'm going to guess is possible by the fact that lots of square roots will cancel, but there's no doubt it will be ugly. This is my solution, and see if you can finish it off and rewrite this ugly mess I got at the end.
Cheers.
EDIT: ONE OF THE BOUNDS IS WRONG ON THE SPHERICAL TRIPLE INTEGRAL. Instead of $0$ to $\phi = 2(\arctan(\frac{\sqrt{2 - b^2 + 2\cos\theta \sin\theta}}{b - 1}))$, you need to use $\phi = 2(\arctan(\frac{\sqrt{2 - b^2 + 2\cos\theta \sin\theta}}{b - 1}))$ to $\frac{\pi}{2}$. Adjust the integral accordingly. Still call this expression $f(\theta)$, but the final integral that I wrote is slightly different. Plug in $f(\theta)$ and $\frac{\pi}{2}$ instead of $0$ and $f(\theta)$.
Edit with proposed solution:
Please double check my work, but this is the solution I found.
Our integral, with the corrected spherical coordinate bounds, and all of the constants put together, becomes:
$$\frac{1}{\sqrt{\pi}} \int_{0}^{\frac{\pi}{4}} f(\theta)d\theta$$
I'm going to rename the variable of integration to $x$, and replace what we decided $f(\theta)$ is, to get the following:
$$\frac{1}{\sqrt{\pi}} \int_{0}^{\frac{\pi}{4}} \mathrm{cos}(2\arctan(\frac{\sqrt{2 - b^2 + 2\cos x \sin x}}{b - 1}))dx.$$
I will leave this for you to verify, but this cosine bit should turn into:
$$\mathrm{cos}(2\arctan(\frac{\sqrt{2 - b^2 + 2\cos x \sin x}}{b - 1})) = \frac{2b^2-2b-1-\sin (2x)}{2b+1+ \sin (2x)}$$
Just using the cosine double angle, and simplifying some fractions. Integrating this from $0$ to $\frac{\pi}{4}$ and then a simple u-substitution, yields:
$$\int_{0}^{\frac{\pi}{4}} \cos (2\arctan(\frac{\sqrt{2 - b^2 + 2\cos x \sin x}}{b - 1})) dx = \int_{0}^{\frac{\pi}{2}} \frac{2b^2-2b-1- \sin (x)}{2b+1+ \sin (x)}dx$$.
We'll use some shorthand, and let $\alpha = 2b^2-2b-1$ and $\beta = 2b+1$. This integral gets messy again, so this will become helpful.
Again, I apologize for this, but the answer to this integral in terms of $\alpha$ and $\beta$ is:
$$-\dfrac{\left(2{\beta}+2{\alpha}\right)\sqrt{1-{\beta}^2}\ln\left(\frac{\left|2\sqrt{1-{\beta}^2}-2{\beta}+2\right|}{\left|2\sqrt{1-{\beta}^2}+2{\beta}-2\right|}\right)+\left(-2{\beta}-2{\alpha}\right)\sqrt{1-{\beta}^2}\ln\left(\frac{2\sqrt{1-{\beta}^2}+2}{\left|2\sqrt{1-{\beta}^2}-2\right|}\right)+{\pi}{\beta}^2-{\pi}}{2\left({\beta}^2-1\right)}$$
Replace all of this with what we let $\alpha$ and $\beta$ equal to, and you'll get something monstrous, or maybe not! Some of these $\alpha + \beta$ expressions and $\beta ^2 - 1$ expressions may reduce to something nicer. The point is, it's possible. And don't forget to add your $\frac{1}{\sqrt{\pi}}$ out the front. I also understand that we've let some more restrictions on $\beta$, and hence $b$, in addition to the restriction we set on $b$ during the set up of the triple integral that I briefly mentioned, but this is not my place to start figuring out how to mitigate that, or whether it matters, or anything like that, but this is the answer you get from me. My head hurts. The point is: it's possible. Is it nice? Aboslutely not.
Cheers, and I hope this answer is helpful. I recommend you do this by hand, as awful as that sounds, to double-check my integration and logic.