[proofplan]
We prove the existence of each [partial derivative](/page/Partial%20Derivative) by restricting the differentiability approximation to the coordinate line through $a$ in the direction $e_i$. Openness of $U$ ensures that this coordinate line remains inside $U$ for all sufficiently small parameters. Dividing the differentiability expansion by the parameter gives the desired directional difference quotient, and the remainder vanishes because $|t e_i|=|t|$. Finally, applying coordinate projections identifies the entries of the Jacobian matrix with the scalar component partial derivatives.
[/proofplan]
[step:Restrict the differentiability expansion to a coordinate line]
Fix $i \in \{1,\ldots,m\}$. Since $U$ is open and $a \in U$, there exists $\rho>0$ such that $B(a,\rho)\subset U$. Hence, for every $t \in \mathbb{R}$ with $0<|t|<\rho$, the point $a+t e_i$ belongs to $U$ because $|t e_i|=|t|<\rho$.
Since $f$ is differentiable at $a$, there is a map $r: B(0,\rho)\to \mathbb{R}^n$ such that, for every $h \in B(0,\rho)$,
\begin{align*}
f(a+h)=f(a)+Df_a(h)+r(h),
\end{align*}
and
\begin{align*}
\lim_{h\to 0}\frac{|r(h)|}{|h|}=0.
\end{align*}
Substituting $h=t e_i$ gives, for every $t \in \mathbb{R}$ with $0<|t|<\rho$,
\begin{align*}
f(a+t e_i)-f(a)=Df_a(t e_i)+r(t e_i).
\end{align*}
[/step]
[step:Take the coordinate difference quotient and pass to the limit]
Because $Df_a:\mathbb{R}^m\to\mathbb{R}^n$ is linear, for every $t\ne 0$,
\begin{align*}
\frac{Df_a(t e_i)}{t}=Df_a(e_i).
\end{align*}
Therefore, for $0<|t|<\rho$,
\begin{align*}
\frac{f(a+t e_i)-f(a)}{t}=Df_a(e_i)+\frac{r(t e_i)}{t}.
\end{align*}
The remainder term tends to $0$ in $\mathbb{R}^n$ because
\begin{align*}
\left|\frac{r(t e_i)}{t}\right|=\frac{|r(t e_i)|}{|t|}=\frac{|r(t e_i)|}{|t e_i|}\to 0
\end{align*}
as $t\to 0$. Hence the coordinate-direction limit exists and equals $Df_a(e_i)$:
\begin{align*}
\partial_{x_i}f(a)=\lim_{t\to 0}\frac{f(a+t e_i)-f(a)}{t}=Df_a(e_i).
\end{align*}
[guided]
The definition of the partial derivative in the $x_i$ direction asks for the limit of the one-variable difference quotient obtained by moving from $a$ only in the coordinate direction $e_i$. The first point to check is that this quotient is defined. Since $U$ is open and $a\in U$, there is a radius $\rho>0$ such that $B(a,\rho)\subset U$. If $0<|t|<\rho$, then $|a+t e_i-a|=|t e_i|=|t|<\rho$, so $a+t e_i\in U$.
Differentiability at $a$ means that the increment of $f$ is the linear increment $Df_a(h)$ plus an error that is small compared with $|h|$. Thus there is a map $r:B(0,\rho)\to\mathbb{R}^n$ such that
\begin{align*}
f(a+h)=f(a)+Df_a(h)+r(h)
\end{align*}
for $h\in B(0,\rho)$, and
\begin{align*}
\lim_{h\to 0}\frac{|r(h)|}{|h|}=0.
\end{align*}
Now choose the special increment $h=t e_i$. This gives
\begin{align*}
f(a+t e_i)-f(a)=Df_a(t e_i)+r(t e_i).
\end{align*}
Dividing by $t$ is legitimate for $t\ne 0$, and linearity of $Df_a$ gives
\begin{align*}
\frac{Df_a(t e_i)}{t}=Df_a(e_i).
\end{align*}
Therefore
\begin{align*}
\frac{f(a+t e_i)-f(a)}{t}=Df_a(e_i)+\frac{r(t e_i)}{t}.
\end{align*}
The remaining term vanishes as $t\to 0$. Indeed, since $|t e_i|=|t|$,
\begin{align*}
\left|\frac{r(t e_i)}{t}\right|=\frac{|r(t e_i)|}{|t|}=\frac{|r(t e_i)|}{|t e_i|}\to 0.
\end{align*}
Thus the coordinate difference quotient converges in $\mathbb{R}^n$, and its limit is exactly $Df_a(e_i)$. By the definition of the vector-valued partial derivative,
\begin{align*}
\partial_{x_i}f(a)=Df_a(e_i).
\end{align*}
[/guided]
[/step]
[step:Identify the component partial derivatives and the Jacobian matrix]
Let $\pi_j:\mathbb{R}^n\to\mathbb{R}$ denote the $j$-th coordinate projection, so that $f_j=\pi_j\circ f$ for each $j\in\{1,\ldots,n\}$. From the vector limit just proved, applying the continuous [linear map](/page/Linear%20Map) $\pi_j$ gives
\begin{align*}
\partial_{x_i}f_j(a)=\pi_j(\partial_{x_i}f(a))=\pi_j(Df_a(e_i)).
\end{align*}
Thus the $i$-th column of the matrix of $Df_a$ in the standard bases is the vector $\partial_{x_i}f(a)$, and its $j$-th entry is $\partial_{x_i}f_j(a)$. Hence
\begin{align*}
(Jf_a)_{ji}=\partial_{x_i}f_j(a),
\end{align*}
so $Jf_a$ is precisely the standard-basis matrix representation of $Df_a$.
[/step]
