When considering a function $f$ from some open region in $\mathbb{R^n}$ some other region in $\mathbb{R^m}$, I have seen the derivative of $f$ at $x$ defined as the $m\times n$ matrix $df$ (a.k.a the Jacobian) such that:
$$df(x) \xi = \lim_{t\to 0} \frac{f(x+t\xi)-f(x)}{dt},\ \forall\xi\in\mathbb{R}^n$$ My question is, how can we prove that this limit exists and is the same for all $\xi$? In other words, suppose we found the derivative for a given $\xi$, how do we know this derivative is the same for other $\xi$ as well?
Pretty sure I was taught this back in the day...
Most of the time one does not use the definition to prove that the derivative exists. Instead, one uses a property such as: $f$ is the sum or the product or a linear combination or a composite or a fraction etc of derivable functions and there are theorems for this cases.
On singular points, the directional derivatives are a way to reduce the problem to the well known question of computing the derivative of a function of only one real variable. The result is a limit depending on $\xi$
$$\ell(\xi) = \lim_{t\to 0}\frac{f(x+t\xi)-f(x)}{t}$$ Then one must study the map $\xi \mapsto\ell(\xi)$ from $\mathbb{R}^n$ into $\mathbb{R}^m$. If this map turns out to be linear, then $f$ is derivable at $x$. There are examples for which the map $\ell$ exists, but it is not linear. For example $$ f(x, y) = \frac{x|y|}{\sqrt{x^2+y^2}} $$ with $f(0,0)=0$ satisfies this property.
As said in the comments, a common situation that guarantees the existence of the derivative is when the partial derivatives exist in a neighbourhood of $x$ and are continuous at $x$. A sketch of the proof goes as follows: let $d^i\in \mathbb{R}^m$ be the $i$-th first-order partial derivative of $f$ at point $x$, and let $e^1,\ldots,e^n$ be the standard basis of $\mathbb{R}^n$. Let $\xi = \sum \xi_k e^k \in \mathbb{R}^n$. We can define by induction a sequence $x^k$ by $x^0=x$ and $x^{k}=x^{k-1}+\xi_k e^k$, so that $x^n=x+\xi$. Finally, let us define $g(y) =f(y) - \sum d^i (y_i-x_i)$, so that $g$ has continuous partial derivatives in a neighbourhood of $x$ and this derivatives vanish at point $y=x$. We then write $$ f(x+\xi) -f(x)-\sum d^i \xi_i = (g(x^n)-g(x^{n-1}))+\ldots + (g(x^1)-g(x^0)) $$ Then one proves by the mean value theorem that $$ |g(x^{k})-g(x^{k-1})|\le |x^{k}-x^{k-1}| \sup_{|z-x|\le |\xi|} |\frac{\partial_k}{\partial x_k} f(z) - d^k | $$ This proves that $$ |f(x+\xi) -f(x)-\sum d^i \xi_i | \le |\xi| \epsilon(\xi) $$ where $\epsilon(\xi)\to 0$ when $\xi\to 0$. It follows easily that $f$ is derivable at $x$ with derivative $df(x)\xi=\sum d^i \xi_i $.