Inspired by the Bürmann's expansion of the error function find on wikipedia I found the begin of an expansion also using $e^{-x^2}$ :
Explanation for $x\in(0,1]$:
We start with a well-know inequality for $x\geq 0$ :
$$\operatorname{erf}\left(x\right)-\tanh\left(x\right)\geq 0$$
The problem is the maximum is not an integer easily reachable so using a limit we divide by $x$ to get :
$$g(x)=\frac{\operatorname{erf}\left(x\right)-\tanh\left(x\right)}{x}$$
The maximum define as a limit is at $x\to 0$
Now we use a famous function $f(x)=e^{-x^2}$ .For that we multiply the function $f$ by the limit a zero of $g$ and we substract this quantity to $g(x)$ we get :
$$h(x)=\frac{\operatorname{erf}\left(x\right)-\tanh\left(x\right)}{x}-\left(-1+\frac{2}{\sqrt{\pi}}\right)e^{-x^{2}}$$
Again the maximum is not easily reachable so now we divide by $x^2$, substract and we have :
$$r\left(x\right)=\frac{\frac{\operatorname{erf}\left(x\right)-\tanh\left(x\right)}{x}-\left(-1+\frac{2}{\sqrt{\pi}}\right)e^{-x^{2}}}{x^{2}}-e^{-x^{2}}\left(-\frac{2}{3}+\frac{\frac{4}{3}}{\sqrt{\pi}}\right)$$
As I have some problem with Desmos I have only this .
In fact I can pursue it the next step is :
$$\frac{\frac{\frac{\operatorname{erf}\left(x\right)-\tanh\left(x\right)}{x}-\left(-1+\frac{2}{\sqrt{\pi}}\right)e^{-x^{2}}}{x^{2}}-e^{-x^{2}}\left(-\frac{2}{3}+\frac{\frac{4}{3}}{\sqrt{\pi}}\right)}{x^{2}}-\left(-\frac{3}{10}+\frac{8}{15}\cdot\frac{1}{\sqrt{\pi}}\right)e^{-x^{2}}$$
How to finish it and justify it ?