These days I've been looking for a rigurous proof of the multivariable chain rule and I've finally found one that I think is very easy to understand. I will leave it here (if nobody minds) for anybody searching for this that is not familiar with little-o notation, Jacobians and stuff like this. To understand this proof, all you need to know is the mean value theorem.
Let's say we have a function $f(x,y)$ and $x = x(t), y = y(t)$. Let's also take $z(t) = f(x(t), y(t))$ By definition, the derivative of z $z'(t)$ is
$$ z'(t) = \lim_{\Delta t \to 0}{\frac {f(x(t+\Delta t),y(t+\Delta t)) - f(x,y)}{\Delta t}}$$.
$$ Let \ \Delta x = x(t+\Delta t)-x(t),$$ $$\Delta y = y(t+\Delta t)-y(t)$$
Now I'll take the numerator of the fraction in the limit, and make a small change.
$$ f(x(t+\Delta t), y(t+\Delta t)) - f(x,y) = f(x+\Delta x, y+\Delta y) - f(x,y)$$
$$ = \left[f(x+\Delta x, y+\Delta y) - f(x+\Delta x, y)\right] + \left[f(x+\Delta x, y) - f(x, y)\right]$$
I have just added and substracted $f(x+\Delta x, y)$. For some reason, I will invert the terms.
$$ = \left[f(x+\Delta x, y) - f(x, y)\right] + \left[f(x+\Delta x, y+\Delta y) - f(x+\Delta x, y)\right]$$.
Now, let's define 2 functions and I will name them g and h. First,
$$ Let \ g(x) = f(x, y) \implies g'(x) = \frac {\partial f} {\partial x} $$.
Please note that y is constant here since g is a function of a single variable. Now, by the mean value theorem we have
$$ \exists c_1 \in (x, x+\Delta x) \ so \ that$$
$$\frac {g(x+\Delta x) - g(x)} {\Delta x} = g'(c_1) $$
$$ \Longleftrightarrow $$
$$ f(x+\Delta x, y) - f(x, y) = f_x(c_1, y)\Delta x$$
Similarly, using the function
$$ h(y) = f(x + \Delta x, y) \implies h'(y) = \frac {\partial} {\partial y}f(x+\Delta x, y)$$
We will have by the same logic that
$$ f(x+\Delta x, y + \Delta y) - f(x+\Delta x, y) =
f_y(x + \Delta x, c_2)\Delta y, c_2 \in (y, y+\Delta y) $$
Notice that $c_1$ and $c_2$ are bounded with respect to $\Delta x$ and $\Delta y$
So as $\Delta x \to 0, c_1 \to x$ and as $\Delta y \to 0, c_2 \to y$. By our definition of $\Delta x$ and $\Delta y$, as $\Delta t \to 0$, both $\Delta x$ and $\Delta y$ $\to 0$. So, as $\Delta t \to 0$, $c_1 \to x$ and $c_2 \to y$.
The last step of the proof is to sum this all up, divide by $\Delta t$ and take the limit as $\Delta t \to 0$
$$ f(x(t+\Delta t), y(t+\Delta t)) - f(x, y) = f_x(c_1, y)\Delta x + f_y(x+\Delta x, c_2)\Delta y $$
$$ \lim_{\Delta t \to 0} \frac {f(x(t+\Delta t), y(t+\Delta t))}{\Delta t} = \lim_{\Delta t \to 0} f_x(c_1, y)\frac {\Delta x}{\Delta t} + f_y(x+\Delta x, c_2)\frac {\Delta y}{\Delta t} = f_x(x, y)x'(t) + f_y(x, y)y'(t) \ QED $$
Edit: After a long time I've realised that this proof assumes that $f$ has partial derivatives defined on intervals around the point $(x, y)$ and they are continuous at the point. This is a sufficient condition for the function to be ($\mathbb{R}^2$-)differentiable at $(x, y)$, but it's not equivalent. Yet, the multivariable chain rule works for the function being just differentiable at that point. So for a general proof, one should first understand little-o notation as in the other answers.
The "obvious" formula to induct (assuming your formula in 2D is right) to me would be that: there exists $c=(c_1\dots c_n)$ such that $$\frac{\partial^nf}{\partial x_1\dots \partial x_n}(c) \ (b_n-a_n)\dots (b_1-a_1) = \delta_1(\dots (\delta_n f)\dots))$$
where $\delta_k g := g(b_k)-g(a_k)$ only changes the inputs in the $k$th variable. Writing $\partial_i$ for $\partial/\partial x_i$ and $V_n(a,b):=(b_n-a_n)\dots (b_1-a_1) $ we have
\begin{align} \partial_1f(c) V_1(a,b) &= \delta_1 f=f(b_1)-f(a_1)\\
\partial_1\partial_2f(c) V_2(a,b) &= \delta_1 \delta_2 f =\delta_1(f(\bullet,b_2)-f(\bullet,a_2)) \\&= f(b_1,b_2)-f(b_1,a_2) - [f(a_1,b_2)-f(a_1,a_2)]\end{align}
matching the two earlier formulas, and for your convenience here's the third one:
\begin{align} \partial_1\partial_2\partial_3f(c) V_3(a,b)&= \delta_1\delta_2\delta_3 f=\delta_1(f(\bullet , b_2,b_3)-f(\bullet , b_2,a_3) - [f(\bullet , a_2,b_3)-f(\bullet , a_2,a_3)])\\ &=
(f(b_1 , b_2,b_3)-f(b_1, b_2,a_3) - [f(b_1 , a_2,b_3)-f(b_1, a_2,a_3)])
\\
&\quad -
(f(a_1, b_2,b_3)-f(a_1 , b_2,a_3) - [f(a_1 , a_2,b_3)-f(a_1 , a_2,a_3)]) \end{align}
Best Answer
$\def \R {\mathbb{R}}$ It depends on what you mean by mean value theorem in several variables. What doesn't work is mean value theorem for $f: \R^n \to \R^m$ for $m > 1$ since each coordinate in codomain can dictate a different point in domain. But the case $m = 1$, $n > 1$ is ok.
The conterexample on page 2 of InvFT is a counterexample for following mean value theorem: For differentiable $f: \R^n \to \R^m$, there is a point $c$ on line $[a, b] ⊆ \R^n$ such that $f(b) - f(a) = (D_c f)(b - a)$, where $D_c f: \R^n \to \R^m$ is the derivative / differential of $f$ at point $c$. However this theorem holds if $m = 1$ as Theorem 9 from your scanned source shows.
But Theorem 11 on page 14 from your scanned source says something different. It says that there are points $c_i$ on line $[a, b]$, $i ∈ \{1, …, m\}$ such that $f(b) - f(a) = L(b - a)$ where $L = [D_{c_i} f_i: i ∈ \{1, …, m\}]$ where $f_i$ is $i$-th component of $f$.