$\int\text{e}^{-ax^2 } \text{erf}\left(bx + c\right) dx$

I'm hoping to find a closed expression for the following integral. $$ \int\text{e}^{-ax^2 } \text{erf}\left(bx + c\right) dx $$ One can find a solution for a family of products between exponentials and error functions. None of which apparently have the offset term in the error function.

I have tried tackling the problem with two approaches.

Approach #1: Expanding the error function hoping to find nice cancelations leading to the maclerin series of some known elementary function. Following a similar approach by Alex:

$$ \begin{aligned} \int\text{e}^{-ax^2 } \text{erf}\left(bx + c\right) dx &= \frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n }{n!(2n+1)} \int(bx+c)^{2n +1} \text{e}^{-ax^2} dx \\ & = \frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n }{n!(2n+1)} \int \sum_{k=0}^{2n+1} {{2n+1}\choose{k}} (bx)^{k} c^{2n+1-k} \text{e}^{-ax^2} dx \\ & = \frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n }{n!(2n+1)} \sum_{k=0}^{2n+1} {{2n+1}\choose{k}} c^{2n+1-k} b^k\int x^{k} \text{e}^{-ax^2} dx = \\ &\frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n }{n!(2n+1)} \sum_{k=0}^{2n+1} {{2n+1}\choose{k}} c^{2n+1-k} b^k \left(-\frac{1}{2}a^{-\frac{k+1}{2}} \Gamma\left(\frac{k+1}{2},ax^2\right)\right) \\ &= -\frac{1}{\sqrt{a}\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n }{n!(2n+1)} \sum_{k=0}^{2n+1} {{2n+1}\choose{k}} c^{2n+1-k} \left(\frac{b}{\sqrt{a}}\right)^k \Gamma\left(\frac{k+1}{2},ax^2\right) \end{aligned} $$

I have used the binomial expansion for $(bx + c)^{2n+1}$ and that $\int x^k \text{e}^{-ax^2}dx = -\frac{1}{2}a^{-\frac{k+1}{2}} \Gamma\left(\frac{k+1}{2} ,ax^2\right)$ where $\Gamma(,)$ is the incomplete gamma function. Too bad, the last term can not be combined again in the form of a binomial expansion.

Approach #2: Instead of expanding the error function, I tried writing it in terms of the cumulative CDF function (Q-Function) as $\text{erf}(x) = 2\Phi(\sqrt{2} x) - 1$. However, the following can be shown to be true using integration under the integral sign with respect to $\mu$. [Section 2.4, and ref] $$ \frac{1}{\sqrt{2 \pi} \sigma}\int_{-\infty}^{\infty}\Phi(\lambda x) \text{e}^{-\frac{(x - \mu)^2}{2 \sigma^2}}dx =\Phi\left(\frac{\lambda \mu}{\sqrt{1+\lambda^2\sigma^2}}\right) $$

Now with some change of variables and rescaling we are instead interested in the following integral: $$ \int\text{e}^{a_1 x^2 + a_2 x} \text{erf}\left(x\right) dx = \underbrace{2\int \text{e}^{a_1 x^2 + a_2 x} \Phi(\sqrt{2} x) dx}_{I} - \underbrace{\int \text{e}^{a_1 x^2 + a_2 x} dx}_{easy} $$

However, what I'm not certain of if I can use the trick of integration under integral sign for the indefinite integral labeled I. Can I, with some change of variables, use the result deduced for the definite integral case as $\Phi\left(\frac{\lambda \mu}{\sqrt{1+\lambda^2\sigma^2}}\right) + C(x)$?

EDIT: It seems that the problem in had has no closed form solution as pointed out by user90369. Also, user90369 has pointed out that the following more general case have no closed form solution. $$ \int x^{2n} \text{e}^{-ax^2} \text{erf}(bx+c) dx $$ I was wondering, if there are any good approximations that I can use here. By good, I mean refer to an error that is $|e(x)| \leq 10^{-5} \forall x$. For starter, I was looking at the high accuracy approximations in here for the erf function. Unfortunately, none of these approximations result into an integral that inherits a closed form solution. I, however, have the following suggested approach with the use of the following identity. $$ \text{erf}(bx+c) = 2 \Phi\left(\sqrt{2} (bx+c)\right) - 1 $$ This results into the following: $$ \begin{aligned} \int\text{e}^{-ax^2 } \text{erf}\left(bx + c\right) dx = 2\int \text{e}^{-ax^2 } \Phi\left(\sqrt{2} (bx+c)\right) dx - \int \text{e}^{-ax^2 } dx \end{aligned} $$ Now, one can use the approximation of the $\Phi$-function that results from applying Chernof's bound. Link $$ \Phi(x) \approx \frac{1}{12} \text{e}^{-\frac{x^2}{2}} + \frac{1}{4} \text{e}^{-\frac{2}{3} x^2} $$ I'd like to take a suggestion of how good is this approximation after computing the integral. Or maybe if there are other better approximations/recommendations that result into a manageable integral afterwards.

I don't know if this helps for an useful approximation but maybe it's better than nothing. :-)

For $\,v\in\mathbb{N}_0\,$ we get

$$ \int x^{2v+1}e^{-ax^2}dx= -\frac{v!e^{-ax^2}}{2a^{v+1}}\sum\limits_{j=0}^v\frac{(ax^2)^j}{j!} + C_{2v+1} $$

and

$$ \int x^{2v}e^{-ax^2}dx= \frac{(2v)!\sqrt{a\pi}\text{erf}(\sqrt{a}x)}{2^v v!(2a)^{v+1}}-e^{-ax^2}\sum\limits_{j=0}^{v-1}\frac{(v-j)!(2v)!x^{2v-2j-1}}{2^j v!(2v-2j)!(2a)^{j+1}} + C_{2v} $$

and it follows:

\begin{align} & \hphantom{ {}={}} \int e^{-ax^2} \text{erf}(bx+c)dx \\ &= \sum\limits_{k=0}^\infty\frac{(-1)^k}{k!(2k+1)}\int (bx+c)^{2k+1} e^{-ax^2} dx \\ &= \sum\limits_{k=0}^\infty \frac{(-1)^k}{k!(2k+1)}\sum_{v=0}^{2k+1}\binom {2k+1} v b^v c^{2k+1-v}\int x^v e^{-ax^2} dx \\ &= \sum\limits_{k=0}^\infty \frac{(-1)^k}{k!(2k+1)}\sum_{v=0}^k\binom {2k+1} {2v} b^{2v} c^{2k+1-2v}\int x^{2v} e^{-ax^2} dx \\ &\hspace{5mm} +\sum\limits_{k=0}^\infty \frac{(-1)^k}{k!(2k+1)}\sum_{v=0}^k\binom {2k+1} {2v+1} b^{2v+1} c^{2k-2v}\int x^{2v+1} e^{-ax^2} dx \\ &= \sum\limits_{k=0}^\infty\frac{(-1)^k}{k!(2k+1)} \sum\limits_{v=0}^k \binom {2k+1} {2v} b^{2v}c^{2k-2v+1} \\ &\hspace{3cm} \cdot \left( \frac{(2v)!\sqrt{a\pi}\text{erf}(\sqrt{a}x)}{2^v v!(2a)^{v+1}}-e^{-ax^2}\sum\limits_{j=0}^{v-1}\frac{(v-j)!(2v)!x^{2v-2j-1}}{2^j v!(2v-2j)!(2a)^{j+1}} \right) \\ &\hspace{5mm} - \sum\limits_{k=0}^\infty\frac{(-1)^k}{k!(2k+1)} \sum\limits_{v=0}^k \binom {2k+1} {2v+1} b^{2v+1}c^{2k-2v} \frac{v!e^{-ax^2}}{2a^{v+1}}\sum\limits_{j=0}^v\frac{(ax^2)^j}{j!} + C \\ &= \sqrt{a\pi}\text{erf}(\sqrt{a}x)\sum\limits_{k=0}^\infty\frac{(-1)^k}{k!(2k+1)}\sum\limits_{v=0}^k \binom {2k+1} {2v} \frac{(2v)!b^{2v}c^{2k-2v+1}}{2^v v!(2a)^{v+1}} \\ &\hspace{5mm} -e^{-ax^2}\sum\limits_{k=0}^\infty\frac{(-1)^k}{k!(2k+1)}\sum\limits_{v=0}^k \binom {2k+1} {2v} b^{2v}c^{2k-2v+1}\sum\limits_{j=0}^{v-1}\frac{(v-j)!(2v)!x^{2v-2j-1}}{2^j v!(2v-2j)!(2a)^{j+1}} \\ &\hspace{5mm} -e^{-ax^2}\sum\limits_{k=0}^\infty\frac{(-1)^k}{k!(2k+1)}\sum\limits_{v=0}^k \binom {2k+1} {2v+1} b^{2v+1}c^{2k-2v} \frac{v!}{2a^{v+1}}\sum\limits_{j=0}^v\frac{(ax^2)^j}{j!} + C \end{align}