The main takeaway here is that the derivative of functions in $\R^n$ can be defined as linear transformations.
Def 1.7.3. Partial derivatives. Let $U$ be an open subset of $\R^n$ and $f:U\to \R$ a function. The partial derivative of $f$ with respect to the $i$th variable, and evaluated at $a$, is the limit
$$ D_if(a) = \lim_{h\to 0}\frac{1}{h}\begin{pmatrix}\begin{pmatrix}a_1\\\vdots\\a_i+h\\\vdots\\a_n\end{pmatrix} - f\begin{pmatrix}a_1\\\vdots\\a_i\\\vdots\\a_n\end{pmatrix}\end{pmatrix}, $$
if the limit exists.
Def 1.7.7. The Jacobian Matrix. Let $U$ be an open subset of $\Reals^n$. The Jacobian Matrix of a function $f: U \mapsto \Reals^m$ is the $m\times n$ matrix composed of the $n$ partial derivatives of $f$ evaluated at $a$:
$$ \begin{bmatrix}\mathbf{Jf(a)}\end{bmatrix} = \begin{bmatrix} D_1f_1(a) & \dots & D_nf_1(a)\\ \vdots & \ddots & \vdots\\ D_1f_m(a) & \dots & D_nf_m(a) \end{bmatrix} $$
Prop and Def 1.7.9. The Derivative. Let $U \subset \R^n$be an open subset and let $f: U \to \R^m$ be a mapping; let $a$ be a point in $U$. If there exists a linear transformation $L: \R^n \to \R^m$ such that
$$ \lim_{\vec{\mathbf{h}} \to \vec{\mathbf{0}}} \frac{1}{|\vec{\mathbf{h}}|} \mathbf{((f(a+\vec{h}) - f(a)) - (}L(\mathbf{\vec{h}})) = \vec{\mathbf{0}}, $$
then $f$ is differentiable at $a$, and $L$ is unique and is the derivative of $f$ at $a$, denoted $\mathbf{[Df(a)]}$.
Theorem. If $f$ is differentiable at $a$, then all partial derivatives of $f$ at $a$ exist, and the matrix representing $[Df(a)]$ is $[Jf(a)]$.
Proof of 1.7.7. We prove that the Jacobian Matrix is a unique linear transformation. Because the linear transformation $L$ is represented by the matrix whose $i$th column is $L(\vec{\mathbf{e_i}})$, we show that
$$ L(\vec{\mathbf{e_i}}) = \vec{D_i}f(a), $$
where $\vec{D_i}f(a)$ is by definition the $i$th column of the Jacobian matrix $\begin{bmatrix}\mathbf{Jf(a)}\end{bmatrix}$. Using the definition of the derivative (Def 1.7.9), we can set $\vec{\mathbf{h}} = t\mathbf{\vec{e_i}}$ and let $t$ tend to 0.
$$ \lim_{t\vec{\mathbf{e_i}} \to \vec{\mathbf{0}}} \frac{1}{|t\vec{\mathbf{e_i}}|} \mathbf{((f(a+\mathnormal{t}\vec{e_i}) - f(a)) - (}L(t\mathbf{\vec{e_i}})) = \vec{\mathbf{0}}, $$