You can use my limest function. It is on the file exchange.
>> fun= @(x) (x.^3 - 1)./(x-1)
function_handle with value:
Now use limest. It even provides an estimate of how well it thinks that limit is known.
[L,errest] = limest(fun,1)
Is that correct? l'hopital would tell me of course. Thus, if I differentiate the numerator and the demoninator, we would have 3^x^2/1. At x==1, that is 3.
The symbolic toolbox would agree, but you don't want to see that.
But we can still use the symbolic TB, without use of limit, just using l'hopital...
As expected, it returns 3 as the desired limit.
You can find limest on the file exchange, here:
LIMEST uses an adaptive, multi-order Richardson extrapolation scheme, modified to provide also an estimate of the uncertainty at the extrapolation point, all of my invention.)