Here, after using some transform I have been working on, I got this:
$$y(x) = \frac{(\sigma x + \mu) \cdot ((\sigma x + \mu) + 6) \cdot (1 - e^{-6 (\sigma x + \mu)})}{\sqrt{2 \pi (\sigma x + \mu)}}$$
Simply put, it effectively becomes really close $erfc(x)/2$ after almost $x=\mu+2\sigma/3$ to a level that is very useful to any application and its calculation cost is very small. But before $x=\mu+2\sigma/3$ its less accurate, as it starts with an error of almost 20% at $x=\mu$ but the rapidly vanishes, so the question is how to make it more accurate for $x<\mu+2\sigma/3$?.
Edit 1: this was the original approximation:
$$1-\exp(-6x)\cdot\frac{x(x+6)}{\sqrt{2\pi x}} = 1-\frac{1}{2}\text{erfc}(x)\\$$
$$\text{erfc}(x) = \frac{2}{\sqrt{\pi}}\int_x^\infty e^{-t^2}dt$$
LHS=RHS when x>≈4.95909030730313
If you want to have good approximations of the error function, in the same spirit as Vedder, you can write $$\text{erf}(x)\sim \tanh \left(\sum _{i=1}^n a_i\, x^{2 i-1}\right)$$ and, using Taylor expansions of both sides around $x=0$, identify the coefficients.
This would give $$a_1=\frac{2}{\sqrt{\pi }}\qquad a_2=\frac{2 (4-\pi )}{3 \pi ^{3/2}} \qquad a_3= \frac{96-40 \pi +3 \pi ^2}{15 \pi ^{5/2}}$$ $$a_4=\frac{5760-3360 \pi +532 \pi ^2-15 \pi ^3}{315 \pi ^{7/2}}$$
Looking at the error, its maximum value is $3.37\times 10^{-5}$ for $x=2.42$.
If you consider the norm for $n=4$ $$\Phi=\int_{-\infty}^{+\infty}\Bigg(\text{erf}(x)- \tanh \left(\sum _{i=1}^4 a_i\, x^{2 i-1}\right)\Bigg)^2\,dx=1.55\times 10^{-9}$$
Notice that $\frac{167}{148}$ and $\frac{11}{107}$ are rational approximations of $a_1$ and $a_2$.
If you look here, you will find other good approximations. Using $P_1(x)$ allows to invert the function at the price of a quadratic equation in $x^2$.