Suppose first that \(f\) is differentiable at \(\va\text{.}\) Then there exists a linear map \(L \colon \R^n \to \R^m\) such that
\begin{equation*}
f(\vx)-f(\va)=L(\vx-\va)+\vr(\vx),
\end{equation*}
where
\begin{equation*}
\lim_{\vx \to \va}\frac{\norm{\vr(\vx)}}{\norm{\vx-\va}}=0.
\end{equation*}
Write
\(L=(\ell_{ij})\) in matrix form and
\(\vr(\vx)=\langle r_1(\vx),\ldots,r_m(\vx)\rangle\text{.}\) For
\(\vx \ne \va\text{,}\) define
\begin{equation*}
M_{ij}(\vx)
:=
\ell_{ij}
+
r_i(\vx)\frac{x_j-a_j}{\norm{\vx-\va}^2},
\end{equation*}
and set
\(M_{ij}(\va):=\ell_{ij}\text{.}\) Let
\(M(\vx)=(M_{ij}(\vx))\text{.}\) Then for each component
\(i\text{,}\)
\begin{align*}
\sum_{j=1}^n M_{ij}(\vx)(x_j-a_j) \amp = \sum_{j=1}^n \ell_{ij}(x_j-a_j)\\
\amp\phantom{{}={}} + r_i(\vx)\sum_{j=1}^n \frac{(x_j-a_j)^2}{\norm{\vx-\va}^2}\\
\amp = \sum_{j=1}^n \ell_{ij}(x_j-a_j)+r_i(\vx).
\end{align*}
\begin{align*}
M(\vx)(\vx-\va) \amp = L(\vx-\va)+\vr(\vx)\\
\amp = f(\vx)-f(\va).
\end{align*}
It remains to show that
\(M\) is continuous at
\(\va\text{.}\) For
\(\vx \ne \va\text{,}\)
\begin{equation*}
|M_{ij}(\vx)-\ell_{ij}|
=
\left|r_i(\vx)\frac{x_j-a_j}{\norm{\vx-\va}^2}\right|
\le
\frac{|r_i(\vx)|}{\norm{\vx-\va}}
\le
\frac{\norm{\vr(\vx)}}{\norm{\vx-\va}} \to 0.
\end{equation*}
Hence each entry
\(M_{ij}(\vx)\) tends to
\(\ell_{ij}\text{,}\) so
\(M(\vx) \to L\) as
\(\vx \to \va\text{.}\) Thus
\(M\) is continuous at
\(\va\text{.}\)
Conversely, suppose
\begin{equation*}
f(\vx)-f(\va)=M(\vx)(\vx-\va),
\end{equation*}
where \(M(\vx)\) is continuous at \(\va\text{.}\) Let \(L:=M(\va)\text{.}\) Then
\begin{equation*}
f(\vx)-f(\va)-L(\vx-\va)
=
(M(\vx)-M(\va))(\vx-\va).
\end{equation*}
\begin{equation*}
\frac{\norm{f(\vx)-f(\va)-L(\vx-\va)}}{\norm{\vx-\va}}
\le
\norm{M(\vx)-M(\va)}_{\mathrm{op}}.
\end{equation*}
Since
\(M(\vx)\to M(\va)\text{,}\) the right-hand side tends to
\(0\text{.}\) Hence
\begin{equation*}
\lim_{\vx \to \va}
\frac{\norm{f(\vx)-f(\va)-L(\vx-\va)}}{\norm{\vx-\va}}=0,
\end{equation*}
so
\(f\) is differentiable at
\(\va\text{.}\) By uniqueness of the derivative,
\(Df(\va)=L=M(\va)\text{.}\)