He encontrado una bonita y útil aproximación de la función de error cuadrado $$ \mathrm{erf}^{2}\!\left(x\right)=1-\exp\!\left(-\frac{\pi^{2}}{8}x^{2}\right)+\varepsilon\!\left(x\right). $$
He comprobado numéricamente que el error máximo está limitado por $\left|\varepsilon\!\left(x\right)\right| < 61\cdot10^{-4}$ pero me preguntaron si esto podría ser demostrado de alguna manera analítica. O al menos si el orden del error podría determinarse de tal manera.
La función de error se define como $$ \mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp\left(-t^{2}\right)\,\mathrm{d}t = 2\Phi\!\left(x\sqrt{2}\right)-1, $$ donde $\Phi\!\left(x\right)$ es la función de distribución acumulativa normal.
Y de forma más general: ¿el control numérico no es suficiente en estos casos?