Wikipedia lists
$$\int\limits_{-\infty}^{\infty} \operatorname{erf}(ax+b) \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( \frac{-(x-\mu)^2}{2 \sigma^2}\right) \, dx \\ = \operatorname{erf}\left(\frac{a\mu+b}{\sqrt{1+2a^2\sigma^2}}\right)$$
How is this obtained? Is there an indefinite solution to this integral?
There's an elegant solution based on @fGDu94's hint using probability theory. (I'm sure I've used this argument in a previous question of which this one might be a duplicate, but I can't find it.) Let's first restate the claimed result with $\Phi$ instead of $\operatorname{erf}$ as$$\int_{\Bbb R}(2\Phi((ax+b)\sqrt{2})-1)\frac{d}{dx}\Phi\left(\frac{x-\mu}{\sigma}\right)dx=2\Phi\Bigg((a\mu+b)\sqrt{\frac{2}{1+2a^2\sigma^2}}\Bigg)-1,$$or equivalently as$$\int_{\Bbb R}\Phi((ax+b)\sqrt{2})\frac{d}{dx}\Phi\left(\frac{x-\mu}{\sigma}\right)dx=\Phi\Bigg((a\mu+b)\sqrt{\frac{2}{1+2a^2\sigma^2}}\Bigg).$$The left-hand side is the mean of $\Phi(Y)$ with $Y:=(aX+b)\sqrt{2}$ for $X\sim N(\mu,\,\sigma^2)$ so that $Y\sim N((a\mu+b)\sqrt{2},\,2a^2\sigma^2)$. But in terms of $Z\sim N(0,\,1)$ independent of $Y$ so that $Y-Z\sim N((a\mu+b)\sqrt{2},\,1+2a^2\sigma^2)$,$$\Phi(Y)=P(Z\le Y)=P(Y-Z\ge 0)=\Phi\Bigg((a\mu+b)\sqrt{\frac{2}{1+2a^2\sigma^2}}\Bigg).$$