Note: We need to assume that $f: \mathbb{R}^n \to \mathbb{R}$ is differentiable at a local extremum $\mathbf{a} \in \mathbb{R}^n$ in order to have
$$\lim_{h \to 0} \frac{f(\mathbf{a} + h\mathbf{u} ) - f(\mathbf{a} )}{h}=\lim_{h \to 0^+} \frac{f(\mathbf{a} + h\mathbf{u} ) - f(\mathbf{a} )}{h} \, ,$$
and
$$\lim_{h \to 0} \frac{f(\mathbf{a} + h\mathbf{u} ) - f(\mathbf{a} )}{h}=\lim_{h \to 0^-} \frac{f(\mathbf{a} + h\mathbf{u} ) - f(\mathbf{a} )}{h}\, .$$
Let $\{\mathbf{e}_1, \ldots, \mathbf{e}_n\}$ be the standard basis of $\mathbb{R}^n$.
Suppose $f: \mathbb{R}^n \to \mathbb{R}$ has a local extremum and is differentiable at $\mathbf{a} \in \mathbb{R}^n$.
WLOG, assume $f$ has a local maximum at $\mathbf{a}$. This means there is $\delta>0$ so that $f(\mathbf{x})\leq f(\mathbf{a})$ whenever $\|\mathbf{x}-\mathbf{a}\|<\delta$. So for $j=1, \ldots, n$ we have
$$f(\mathbf{a}+h\mathbf{e}_j)-f(\mathbf{a})\leq 0 \text{ whenever } \|(\mathbf{a}+h\mathbf{e}_j)-\mathbf{a}\|=|h|<\delta.$$
Thus if $h \in (0, \delta)$, then for $j=1, \ldots, n$ we have
$$\frac{f(\mathbf{a}+h\mathbf{e}_j)-f(\mathbf{a})}{h}\leq 0$$
which means $\lim_{h \to 0^+} \frac{f(\mathbf{a} + h\mathbf{e}_j ) - f(\mathbf{a} )}{h} \leq 0. $
Since $\frac{\partial }{\partial x_j}f(\mathbf{a})$ exists $(j=1, \ldots, n)$, we have
$$\frac{\partial }{\partial x_j}f(\mathbf{a})=\lim_{h \to 0^+} \frac{f(\mathbf{a} + h\mathbf{e}_j ) - f(\mathbf{a} )}{h} \leq 0.$$
On the other hand, if $h \in (-\delta,0)$, then for $j=1, \ldots, n$ we have
$$\frac{f(\mathbf{a}+h\mathbf{e}_j)-f(\mathbf{a})}{h}\geq 0$$
which means $\lim_{h \to 0^-} \frac{f(\mathbf{a} + h\mathbf{e}_j ) - f(\mathbf{a} )}{h} \geq 0. $
Since $\frac{\partial }{\partial x_j}f(\mathbf{a})$ exists $(j=1, \ldots, n)$, we have
$$\frac{\partial }{\partial x_j}f(\mathbf{a})=\lim_{h \to 0^-} \frac{f(\mathbf{a} + h\mathbf{e}_j ) - f(\mathbf{a} )}{h} \geq 0.$$
We combine to conclude that $\frac{\partial }{\partial x_j}f(\mathbf{a})=0$, for each of $j=1, \ldots, n$.