With the substitution $x=\frac12(1+\tanh t)$, we obtain $$\int^1_0 e^{-\frac1{c\,x\,(1-x)}}\,dx=\frac12\,\int^\infty_{-\infty}\frac1{\cosh^2 t}\,e^{-\frac4c\,\cosh^2 t}\,dt=\int^\infty_0\frac1{\cosh^2 t}\,e^{-\frac4c\,\cosh^2 t}\,dt,$$ so let's investigate
$$I(x)=\int^\infty_0\frac1{\cosh^2 t}\,e^{-x\,\cosh^2 t}\,dt.$$
Since $$\frac{d}{dt}\frac1{\cosh t}e^{-x\,\cosh^2t}=\left(-\frac{\sinh t}{\cosh^2t}-2x\sinh t\right)e^{-x\,\cosh^2t},$$ we have
$$\int^\infty_0\left(-\frac{\sinh^2 t}{\cosh^2t}-2x\sinh^2 t\right)e^{-x\,\cosh^2t}\,dt=\int^\infty_0\sinh t\,\frac{d}{dt}\frac1{\cosh t}e^{-x\,\cosh^2t}\,dt=-\int^\infty_0e^{-x\,\cosh^2t}\,dt$$ by partial integration. Now we know that
$$1-\frac{\sinh^2 t}{\cosh^2t}=\frac1{\cosh^2t},$$ this means
$$\int^\infty_0\frac1{\cosh^2 t}\,e^{-x\,\cosh^2 t}\,dt=2x\int^\infty_0\sinh^2t\,e^{-x\,\cosh^2 t}\,dt.$$ We make use of the well-known identities
$$\sinh^2t=\frac12(\cosh2t-1)$$ and $$\cosh^2t=\frac12(\cosh2t+1)$$ to get
$$I(x)=x\,e^{-x/2}\int^\infty_0(\cosh2t-1)\,e^{-x/2\,\cosh2t}\,dt.$$ Substituting $t\to t/2$ and taking into account the integral representation $$K_\nu(z)=\int^\infty_0e^{-z\cosh t}\,\cosh(\nu t)\,dt$$ (this is formula 10.32.9 in https://dlmf.nist.gov/10.32), we finally arrive at
$$I(x)=\frac{x}2\,e^{-x/2}(K_1(x/2)-K_0(x/2)).$$
As this almost looked too good to be true, I checked it numerically, with $x=4/c$ and various values of $c$ between $0.1$ and $10$, and found perfect agreement.
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.
Best Answer
Surprisingly the requested integral is easier than the one evaluated. Set $\alpha=z-x~,~ \beta= z-y$. Then with a trivial rescaling of the variable, the integral can be written in the form
$$I(x,y,z)=\frac{1}{T}\int_{0}^1 \frac{du}{{u}^{1/2}(1-u)^{3/2}}\exp\left(-\frac{(\alpha(1-u)+\beta u)^2}{2Tu(1-u)}\right)$$
Expanding the numerator in the exponential and with one more change of variables $q=\frac{|\beta|u}{|\alpha|(1-u)}$ we end up with the form
$$I(x,y,z)=\frac{e^{-\alpha\beta/T}}{T}\sqrt{\frac{|\alpha|}{|\beta|}}\int_{0}^\infty\frac{dq}{\sqrt{q}}\exp\left(- \frac{|\alpha\beta|}{2T}(q+1/q)\right)$$
which can be easily solved for the final result
$$I=\sqrt{\frac{2\pi}{T(z-y)^2}} \exp\left(-\frac{(z-x)(z-y)+|z-x||z-y|}{T}\right)$$