Comment (but probably not solution) too long for the comment field: The series of the exponential gives an alternating series for the given function
$$e^{-x^2}=1-x^2+\tfrac12 x^2-\tfrac16x^4+\tfrac1{24}x^8+\dots+\tfrac{1}{k!}(-x^2)^k+\dots$$
By the quantitative statements of the Leibniz criterion,
$$e_{2n-1}(-x^2)\le e_{2n+1}(-x^2)\le e^{-x^2}\le e_{2n}(-x^2),$$
where $e_n$ is the $n$-th partial sum of the exponential series and the estimate is true for $x^2<2n$. This gives an absolute error of
$$\frac1{(2n)!}x^{4n}\left(1-\tfrac{x^2}{2n+1}\right)\le\left|e_{2n-1}(-x^2)- e^{-x^2}\right|\le\frac1{(2n)!}x^{4n}$$
which shows that in terms of polynomial approximation, one can not do much better than the partials of the exponential series. To get a small relative error, one now needs to chose $n$ large enough to get
$$\frac{e^{15^2}(15)^{4n}}{(2n)!}<ϵ$$
or about
$$2n(\ln(2n)-\ln(225e))>225+|\ln(ϵ)|$$
Which means that you need $n=406,...,419$ to get $ϵ=10^{-4},...,10^{-18}$
Using arithmetic rules of the exponential, one notices that $e^{-x^2}=\left(e^{-(\tfrac{x}{m})^2}\right)^{m^2}$, for instance with $m=15$, and the relative error $ϵ$ is guaranteed if the partial sum $e_n$ is taken such that the relative error of $e_n(-x^2)$ to $e^{-x^2}$ on $[-\tfrac{15}m,\tfrac{15}m]$ is smaller than $\tfrac{ϵ}{2m^2}$. With $m=15$ this requires the easier to control inequality
$$\frac{e^1}{(2n)!}<\tfrac{ϵ}{450}$$
which gives the usual precisions for $n=6,...,12$.
$e_{811}(-15^2)$ on Magma has a gross error, instead of ...e-98 it shows 8e73.
$e_{21}(-\tfrac{x^2}{225})^{225}$ has relative error 5e-19 at $x=15$.
If only basic operations are allowed, then we go back to basic divide and conquer or half-and-square in this case. Take in the above calculus $m$ a nice power of $2$, $m=16$ could work, but let's take $m=32$. Then the error estimate is still largest at the interval end, and the condition at the even larger $x=16$ is
$$\frac{e^{\frac14} \left(\frac14\right)^{2n}}{(2n)!}<\tfrac{ϵ}{2048}
\iff
\frac{\sqrt[4]e}{2^{4n-11}(2n)!}<ϵ$$
Trying out some small values leads to $n=5$ as a likely candidate, and indeed, $e_9(-\tfrac{15^2}{1024})\cdot\exp(15^2)-1\approx -1e-10$
$n=4$ with $e_7$ works as well with an actual error of about $2e-7$, with error bound of actually $1e-6$.
With less than 30 operations, replace 7 by 9 for higher accuracy:
xx=-x*x/1024;
a=xx;
e=1+a;
for k=2 to 7 do a*=xx/k; e+=a; end for;
for k=1 to 10 do e=e*e; end for;
return e
or with a division
xx=x*x/1024;
a=xx;
e=1+a;
for k=2 to 7 do a*=xx/k; e+=a; end for;
e=1/e;
for k=1 to 10 do e=e*e; end for;
return e
First, note that a 4th degree polynomial might not exist, for a simple example, consider finding a quadratic which has points $f(0) = 0, f(1) = 1, f'( \frac{1}{2} ) = 1 $. The first 2 conditions naturally imply that $ f'( \frac{1}{2}) = 0 $, which contradicts the third.
But otherwise, we can use a similar idea to the building blocks of Lagrange Interpolation.
Use $f_1 = B(x-x_2)(x-x_3)(x-x_4)(x-A)$, where $A$ is chosen such that $ f''' (x_5) = 0$ and $B$ is chosen such that $ f_1 (x_1) = 1$. There are no issues here. Define $f_2, f_3, f_4$ similarly.
Use $ f_5 = C (x-x_1)(x-x_2)(x-x_3)(x-x_4) $, where $C$ is chosen such that $ f'''(x_5) = 1$. The only issue here, is if $ [ (x-x_1)(x-x_2)(x-x_3)(x-x_4)]''' (x_5) = 0$, in which case such a constant does not exist. In this case, we need to tag on another linear term, to make a polynomial of degree 6.
Best Answer
The formula that you give for the absolute value of the error is exact, sort of. By that I mean that for every $x$, there is a $\xi(x)$ in our interval such that the error when you use the interpolating quadratic at $x$ is exactly the one obtained by putting $\xi=\xi(x)$ in the error formula.
In most cases, this kind of Mean Value Theorem information leads to estimates, but not to exact values, since we don't know $\xi(x)$.
However, here we are evaluating the third derivative of our cubic at $\xi(x)$. The third derivative is constant! So it doesn't matter what $\xi(x)$ is, that part of the error estimate does not change. The rest is an explicit polynomial in $x$. It follows that the error estimate and the actual error are equal, as you noticed computationally. Both are polynomials in $x$. Given any two polynomials, equality at all values of $x$ (or even at infinitely many values of $x$, or, for cubics, at $4$ values of $x$) means that the coefficients match.
Precisely the same thing happens when you use the same process to approximate a polynomial of degree $n+1$ by a Newton interpolating polynomial of degree $n$. If you look up the error estimate, you will see that it is a certain polynomial in $x$ times a term $f^{(n+1)}(\xi(x))$. (I like that version better than the absolute value version, since we also get information about the sign of the error.)
For a polynomial $f$ of degree $n+1$, the $(n+1)$-th derivative is constant, so
the term $f^{(n+1)}(\xi(x))$ is a constant independent of $x$. Thus the error estimate and the error are identically equal. Since they are both polynomials, that means their coefficients are equal.