I found a nice and useful approximation of squared error function $$ \mathrm{erf}^{2}\!\left(x\right)=1-\exp\!\left(-\frac{\pi^{2}}{8}x^{2}\right)+\varepsilon\!\left(x\right). $$
I checked numerically that maximum error is bounded by $\left|\varepsilon\!\left(x\right)\right| < 61\cdot10^{-4}$ but I was asked if this could be somehow proved analytically. Or at least if order of the error could be determined in such manner.
Error function is defined as $$ \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, $$ where $\Phi\!\left(x\right)$ is normal cumulative distribution function.
And more generally: is numerical check not enough in such cases?