Find an algorithm to calculate the following function I'm solving questions from an exam I failed, and I would love some help with the following question:
Question
We want to calculate the following function in Matlab:
$$
f(x) = \frac{e^{x^2} - (1 + x^2)}{x^4}
$$
We know that as $x$ approaches $0$, $f(x)$ approaches $0.5$, but in Matlab we get: $f(0.0001) = 0$.


*

*Why don't we get a result close to 0.5 for small $x$?

*Propose an algorithm to calculate $f(x)$ for $x = 0.001$ that will be exact for at least 10 decimal points.


My thoughts
I played with the function in Matlab and I saw that when I calculate $1 + x^2$ and $e^{x^2}$ I get 1, which explains the result: x^2 is smaller than the machine's epsilon, so $1+x^2$ cannot be represented in floating point using the number of bits we have, and since both expressions are rounded to 1 we get 0.
I am not sure how to approach 2. I tried to change the form of the expression analytically to something that won't cause similar problems but didn't find anything promising.
Thanks! :)
 A: Let's look how is $e^{x^2} - (1 + x^2)$ evaluated.
First, $e^{x^2}$ is evaluated, and that results in
$$
e^{x^2} = \sum_{k=0}^\infty \frac{x^{2k}}{k!} 
= 1 + x^2 + \frac{x^4}{2} + \dots \approx
1 + 10^{-8} + 0.5 \cdot 10^{-16} + \dots.
$$
But from machine point of view when using double precision arithmetics, $0.5 \cdot 10^{-16}$ is less than unity roundoff, which is $2^{-53} \approx 1.1 \cdot 10^{-16}$. Thus the best double-precision representation for $e^{x^2}$ would be just (nearest 53 bit binary to)
$$
1 + 10^{-8}.
$$
and that is represented precisely the same as $1 + x^2$. Thus the numerator is exactly zero when working with double precision numbers. Dividing zero with $x^4$ effectively does nothing.
For the part 2 i recommend consider
$$
f(x) = \frac{e^{x^2} - (1 + x^2)}{x^4}
$$
as generic power series in $x^2$:
$$
f(x) = \frac{1 + x^2 + \frac{x^4}{2} + \dots - 1 - x^2}{x^4} =
\frac{1}{2} + \frac{x^2}{6} + \dots = \sum_{k=0}^\infty \frac{x^{2k}}{(k+2)!}
$$
I left to you determining of the correct number of the terms needed to represent $f(x)$ with at least 10 digits.
A: Use the Taylor expansion of exponential to get
$$f(x) = \frac{1}{2!} + \frac{1}{3!}x^2 + \frac{1}{4!}x^4 +\frac{1}{5!}x^6 + \dots$$
Matlab gives you $0$ because
$$e^{0.00001}-(1+0.00001) = 0.00000000005000...$$
and maybe it approximates this term to $0$.
