Integrating $\frac1{2\sigma\sqrt{2\pi}}\int_{-\infty}^\infty(1+erf(\frac{z-\mu_i}{\sigma_i\sqrt2}))\exp( -\frac12(\frac{z-\mu}\sigma)^2)dz$

Integral of one error function and Gaussian

Mathematica does not evaluate the integral. Surprisingly, there is an exact result. By changing variables, $z-\mu=x$, your integral of $\operatorname{erf} \times \text{Gaussian}$ may be written as

$$ I(b)=\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(ax+b) \exp \left(-cx^2 \right) $$

Since $\partial_t \operatorname{erf}(t)= \frac{2}{\sqrt{\pi}}e^{-t^2}$, differentiating with respect to $b$ produces

$$ \partial_b I(b)=\int\limits_{-\infty}^\infty dx \ \frac{2}{\sqrt{\pi}} \exp\left( -(ax+b)^2-cx^2 \right) $$

This is just a Gaussian integral over $x$, the result is

$$ \partial_b I(b)= \frac{2}{\sqrt{a^2+c}} \exp \left( -\frac{b^2 c}{a^2 +c} \right) $$

Now we integrate to recover $I(b)$

$$ I(b)-I(0)=\int\limits_0^b \ dt \frac{2}{\sqrt{a^2+c}} \exp \left( - t^2\frac{ c}{a^2 +c} \right) $$

$I(0)=\int dx \ ( \text{odd function})=0$, and the integral on the right is simply the definition of $\operatorname{erf}$, so we have

$$ I(b)=\sqrt{\frac{\pi}{c}} \operatorname{erf} \left( b \sqrt{\frac{c}{a^2+c}} \right) $$

$$ \tag*{$\blacksquare$} $$

Integral of two error functions and Gaussian

We may use a similar method to find that

$$ \int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a_1 x) \operatorname{erf}(a_2 x) e^{-cx^2} = \frac{2}{\sqrt{c\pi}} \tan^{-1}\left( \frac{a_1 a_2}{ \sqrt{c(a_1^2+a_2^2+c)}} \right) $$

However, once we introduce $b\neq 0$, we meet the more complicated integral in an intermediate step

$$ \int\limits_0^b dt \ \operatorname{erf}(\alpha t + \beta) e^{-\gamma t^2} $$

Update! Using the result from this post, we can make progress. Let

$$ J(\beta)=\int\limits_{-\infty}^\infty dx \ \operatorname{erfc}(\alpha x+\beta) \operatorname{erf}(a x+b) \exp \left(-cx^2 \right) $$

Then $J(\infty)=0$. Differentiating

$$ \partial_{\beta} J(\beta)=-\frac{2}{\sqrt{\pi}} \int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a x+b) \exp \left(-cx^2 - (\alpha x +\beta)^2 \right) $$

Which we can evaluate using our first result

$$ \partial_{\beta} J(\beta)=-\frac{2}{\sqrt{\pi}} \sqrt{\frac{\pi}{c+\alpha^2}} \operatorname{erf}\left( \left(b-\frac{a \alpha \beta}{c+\alpha^2} \right) \sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}} \right) \exp \left( -\beta^2 \left( 1+\frac{\alpha^2}{c+\alpha^2}\right) \right) $$

With $\Lambda=\frac{2}{\sqrt{\pi}} \sqrt{\frac{\pi}{c+\alpha^2}}$, $A=-\frac{a \alpha }{c+\alpha^2} \sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}}$, $B=b\sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}}$, and $C=\left( 1-\frac{\alpha^2}{c+\alpha^2}\right)$. Now we integrate

$$ J(\infty)-J(\beta)=-\int\limits_\beta^\infty \ dt \Lambda \operatorname{erf}(At+B)e^{-Ct^2} $$

$$ J(\beta)=\int\limits_\beta^\infty dt \ \Lambda \operatorname{erf}(At+B)e^{-Ct^2}= \frac{\Lambda}{\sqrt{2C}} \int\limits_{\sqrt{2C}\beta}^\infty du \ \operatorname{erf}(Au/\sqrt{2C}+B)e^{-u^2/2} $$

Which may be written in terms of the 'Generalized Owen-T' function

$$ J(\beta)=\frac{2 \sqrt{\pi} \Lambda}{\sqrt{C}}T\left( \sqrt{2C}\beta, A/\sqrt{C},\sqrt{2}B\right) $$

Using identity ii we have

$$ J(\beta)=\frac{2 \sqrt{\pi} \Lambda}{\sqrt{C}} \left[ \frac{1}{2\pi} \left( \tan^{-1}(A/\sqrt{C}) -\tan^{-1}(A/\sqrt{C} + B/\beta\sqrt{C}) \\ -\tan^{-1}(B^{-1}(\beta \sqrt{C} +AB/\sqrt{C} + A^2 \beta /\sqrt{C} )) \right) +\frac{1}{4} \operatorname{erf}(B (1+A^2/C)^{-1/2}) \\ +T(\beta \sqrt{2C},(A\beta^{-1}+B)/\beta \sqrt{C})+T(\sqrt{2}B (1+A^2/C)^{-1/2}, B^{-1}(\beta \sqrt{C} +AB/\sqrt{C} + A^2 \beta /\sqrt{C} ) \right] $$

Where $T$ is the Owen T function. I've checked it numerically, I hope there are no typos here.

Approximate integral of two or more error functions with Gaussian

We can approximate the integral using Laplace's method. I'll demonstrate with two error functions, but in principle you could have more. Let

$$ I(c)=\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a_1 x+ b_1) \operatorname{erf}(a_2 x+ b_2) e^{-cx^2} $$

Expand the product of $\operatorname{erf}$s around the Gaussian peak: $x=0$

$$ \operatorname{erf}(a_1 x+ b_1) \operatorname{erf}(a_2 x+ b_2) = \operatorname{erf}( b_1) \operatorname{erf}( b_2) + x \left( \text{junk} \right) + x^2 \frac{2 e^{-b_1^2-b_2^2}}{\pi} \left( 2 a_1 a_2 -a_2^2 b_2 e^{b_1^2} \sqrt{\pi} \operatorname{erf}(b_1) - a_1^2b_1 e^{b_2^2} \sqrt{\pi} \operatorname{erf}(b_2) \right) + \mathcal{O}(x^3) $$

Term by term integration produces the asymptotic series. The $(\text{junk})$ and all $x^{\text{odd}}$ terms integrate to zero. Noting that

$$ \int\limits_{-\infty}^\infty dx \ x^{2n} e^{-cx^2} = c^{-n-1/2} \Gamma(n+1/2) $$

For fixed $a_i, \ b_i$ we have

$$ I(c) \sim \sqrt{\frac{\pi}{c}} \operatorname{erf}( b_1) \operatorname{erf}( b_2) \\ +\frac{1}{c^{3/2}} \frac{ e^{-b_1^2-b_2^2}}{\sqrt{\pi}} \left( 2 a_1 a_2 -a_2^2 b_2 e^{b_1^2} \sqrt{\pi} \operatorname{erf}(b_1) - a_1^2b_1 e^{b_2^2} \sqrt{\pi} \operatorname{erf}(b_2) \right) \ \ , \ \ c \to \infty $$

Which is the two term approximation.