Interpolation of Gaussian function - minimize relative error I have been trying to interpolate the function $e^{-x^2}$ on interval [-15,15] using standard methods like Lagrange or Newton interpolation for over a month. The goal is for it to be bound by $-\epsilon$ and $\epsilon$. 
No result reached. So I switched to other type of nodes distribution : tried Chebyshev nodes, skewed Chebyshev nodes (see e.g. this link), cosinus nodes... No result reached.
Maybe the skew coefficient was not high enough? I tried to increase the power to which I raise the initial nodes, so they were skewed to the edges even more. No result with different parameters.
Aware of the fact that the higher number of nodes (well, I used about 80-100 nodes in the previous attempts) doesn't necessary lead to higher precision ( because of the Runge's phenomenon), I decided to split the intervals which I work on, and interpolate each of them separately. The edge points were included in the nodes list, because I needed the interpolation polynomial to remain continuous. With the help of this method, I managed to reach the desired error, but absolute, not relative, and also was forced to increase the number of nodes (80 - 45 - 80). 
And right now I am stuck. Could anyone please help, provide an explanation of what I am doing wrong or give me links on any other algorithms (maybe a good explanation of Hermite interpolation or so), so that I am able to implement those methods in Maple? 
P.S. Important note : I am allowed to use only cubic spline. Of course the desired result was reached easily with the spline degree of 4 or 5.
 A: 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

