I had to try my LIMEST tool out on it. As with any adaptive tool, it can be fooled, but it is usually pretty good.
fun = @(x) (exp(x) - 1)./x;
As you can see, fun has problems at zero.
fun(0)
ans =
NaN
although if we evaluate fun near zero, we see it is near 1.
format long g
fun(1e-5)
ans =
1.00000500000696
LIMEST succeeds, and even is able to provide an error estimate of the limit.
[lim,err] = limest(fun,0,'methodorder',3)
lim =
1
err =
2.50668568491927e-15
LIMEST uses a sequence of polynomial approximations, coupled with an adaptive Richardson extrapolation to generate both a limit estimate and a measure of the uncertainty on that limit.
So what problem are you seeing? The failure you saw is simple subtractive cancellation error. Thus, look at the value of
exp(1e-20)
ans =
1
Even with format long g, the double precision value of exp(1e-20) is simply too close to 1 that when we subtract off 1, the result is an exact zero. Divide that by any non-zero value, and we get zero. Of course, when x is actually zero then we have a 0/0 condition, so NaN resulted when I tried that.
Lets see what happens in high precision. I'll use my HPF tool for that computation, and work in 64 decimal digits.
DefaultNumberOfDigits 64
exp(hpf('1e-20'))
ans =
1.000000000000000000010000000000000000000050000000000000000000166
See that when we sutract off 1, the difference between 1 and the exponential value is less than eps(1), so MATLAB must produce a zero value.
exp(hpf('1e-20')) - 1
ans =
1.000000000000000000005000000000000000000016666666666670000000000e-20
The question unasked is how I would choose to generate that function accurately in MATLAB. Clearly, you don't want to use brute force and define fun as I have, as you lose a great deal of accuracy. I'd probably just expand the Taylor series in a limited region around zero, and use fun as above for x significantly different from zero.