There is a small mistake in that you need to use a unit vector. But, other than that, you are right and Wolfram is wrong. What happens is that the computer is using the gradient to calculate the directional derivative. But the formula to calculate the directional derivative using the gradient uses the chain rule, which assumes the function to be differentiable. Your function is not differentiable at $(0,0)$, so the gradient cannot be used to find the directional derivatives there.
It is worth mentioning that WA is making a second mistake: since the partial derivatives are not continuous at $(0,0)$, they cannot be calculated using differentiation rules (and I don't really know how they get the value zero). Calculating explicitly,
$$
\frac{\partial f}{\partial x}f(0,0)=\lim_{h\to 0}\frac{h^3}{h(h^2+0)}=1,
$$
$$
\frac{\partial f}{\partial y}f(0,0)=\lim_{h\to 0}\frac{0}{h(0+h^2)}=0.
$$
Still, the formula $(a,b)\cdot\nabla f(0,0)=a$ gives the wrong result because $f$ is not differentiable.
Some more observations on Mathematica's behavior: it yields
$$ \frac{d}{dx} \int_1^{x} \frac{e^{-t}}{t} \, dt = \frac{e^{-x} - 1}{x} $$
restricted to $\Im(x) \neq 0$ or $\Re(x) \geq 0$. But changing things slightly gives
$$ \frac{d}{dt} \int \frac{e^{-t}}{t} \, dt = \frac{e^{-t}}{t}. $$
I had originally suspected there was something fishy with the branch cut: Mathematica computes
$$ \int_1^x \frac{e^{-t}}{t} \, dt = -\mathrm{Ei}(-1) - \log(x) - \Gamma(0,x) $$
again restricted to $\Im(x) \neq 0 \vee \Re(x) \geq 0$. However:
- The point we are interested in is away from the branch discontinuity
- I would have expected it to get the derivative right even if there were weird branch issues
(using $x+1$ in the above instead of $x$ does not make any qualitative difference)
Without limits, it computes
$$ \int \frac{e^{-t}}{t} \, dt = \mathrm{Ei}(-t) \color{gray}{+ \mathrm{constant}}$$
If you shift the integrand, you get
$$ \frac{d}{dx} \int_0^x \frac{e^{-(u+1)}}{u+1} \, du = \frac{e^{-(u+1)}}{u+1} $$
and correspondingly
$$ \lim_{y \to 0} \frac{ \int_0^y \frac{e^{-(u+1)}}{u+1} \, du }{3 y} = \frac{1}{3e} $$
(I substituted $x^2 = y$ so that wolfram would finish the calculation for me. This substitution does not make a qualitative difference in the original)
I think the key difference is that in the first version, the branch point is $t=0$, and Mathematica focuses on the behavior there -- which is inherently weird and strange because it's a branch point (and given that, I'm not sure if using the result leads to computing something correct but strange, or something ill-defined). But in the second version, the branch point is at $u=-1$, but Mathematica still focuses on $u=0$ so it gets sane results.
Best Answer
WolframAlpha understands the expression $\sqrt[3]{x}$ for negative x in a different way than you expect.
Try this:
lim\frac{\sqrt{1-x}-3}{2+surd(x,3)} as x to -8