# How to perform this Gaussian integral / sum of series?

I am trying to calculate the indefinite integral:

$$\displaystyle\int dx \exp \left(-\dfrac{(a x + b)^2}{2}\right) {\rm erf} \left(\dfrac{x}{\sqrt{2}} \right), \quad x \in \mathbb{R}$$

Context. This integral arises in the integration of a multivariate Normal distribution with finite integration extremes. Let be $$I(x)$$ the integral, so far I have been able to show that, if $$b=0$$:

$$I(x) = \dfrac{1}{a}\sqrt{\dfrac{\pi}{2}} {\rm erf} \left( \dfrac{x}{\sqrt{2}} \right) {\rm erf}\left( \dfrac{a x}{\sqrt{2}} \right) + \\ \quad \dfrac{2}{a \sqrt{\pi}} \exp\left(-\dfrac{x^2}{2} \right) \displaystyle\sum_{k=0}^{\infty} \dfrac{(-1)^k}{k!(2k+1)} \left(\dfrac{a}{\sqrt{2}}\right)^{2k+1} \displaystyle\sum_{j=0}^{k}\dfrac{2k!!}{2j!!}x^{2j}$$

I have done this by expanding the exponential in the integrand in series, then calculating the first four of the resulting integrals, then figuring out the recurrence formulas in k and j. Part of the series can be summed to $${\rm erf}$$, but the second term has proven hard to simplify. Even though the formula above looks like a step in the right direction (to me), I have no idea how to even reach this point if $$b\neq0$$.

For $$b=0$$, the integral can be evaluated exactly in terms of Owen's T function (https://en.wikipedia.org/wiki/Owen%27s_T_function), $$$$I = -\frac{\sqrt{\pi}}{\sqrt{2}a} \left[1+4 T\left(ax,\frac{1}{a}\right) \right].$$$$ For $$b\ne0$$, no idea...