Can we write $f\in C^{1}(\mathbb R^{2})$ as $f(z_{1})-f(z_{2})= (z_{1}-z_{2})\cdot G(z_{1}, z_{2})$? Mean-value theorem for one variable, tells us that if $f:\mathbb R \to \mathbb R$ is continuously differentiable, then we can write, $f(x)-f(y) = (x-y) G(x,y)$; where $x,y \in \mathbb R$ and actually $G(x,y)$ is a constant function.
Now, let $f:\mathbb R^{2} \to \mathbb R^{2}$ such that $f$ is infinitely many times differentiable function on $\mathbb R^{2}.$
For convention, put,  $z_{1}= (x_{1}, y_{1}), z_{2}= (x_{2}, y_{2}) \in \mathbb R^{2}.$

My Question: Can we expect,  $f(z_{1})-f(z_{2})= (z_{1}-z_{2})\cdot G(z_{1}, z_{2})$ ; where $G$ is some function of $z_{1}$ and $z_{2}$ ? If yes, what can we say about  $G$ ?

Thanks,
 A: The following multivariable version of the mean value theorem seems to suit your needs:
Theorem. Assume $f:\>\Omega\to{\mathbb R}^m$ where $\Omega\subset{\mathbb R}^n$ is open, and let $\tilde\Omega:=\{(x,y)\>|\> [x,y]\subset\Omega\}$ be the set of all point pairs $(x,y)$ for which the connecting segment $[x,y]$ lies completely in $\Omega$. Then there is a matrix-valued function 
$$L:\>\tilde\Omega\to R^{m\times n}, \quad (x,y)\mapsto L(x,y)\ ,$$ with
$L(x,x)=df(x)$ (the Jacobian of $f$ at $x$), such that
$$f(y)-f(x)=L(x,y).(y-x)\qquad\forall (x,y)\in\tilde\Omega\ .$$
(The $\ .\ $ denotes application of a linear map, resp. a matrix, to a vector.)
Proof. Consider an $(x,y)\in\tilde\Omega$ and set up the auxiliary function
$$\phi(t):=f\bigl((1-t)x+ty\bigr)\ ,$$
which computes $f$ along the segment $[x,y]$. According to the chain rule one has
$$\phi'(t)=df\bigl((1-t)x+ty\bigr).(y-x)\ .$$ It follows that
$$f(y)-f(x)=\phi(1)-\phi(0)=\int_0^1\phi'(t)\ dt=L(x,y).(y-x)$$
with
$$L(x,y):=\int_0^1 df\bigl((1-t)x+ty\bigr)\ dt\ ,$$
the integration being elementwise on the matrix elements of $df\bigl((1-t)x+ty\bigr)$. The $L$ so constructed is automatically smooth.
