Suppose that the function $f:\mathbb{R}^2\to\mathbb{R}$ has continuous 2nd-order partial derivatives and let $(x_0,y_0)$ be a point in $\mathbb{R}^2$. Prove that for each $\epsilon > 0$ there exists a $\delta > 0$ such that if $0 < |h| < \delta$ and $0 < |k| < \delta$, then $$\left| \frac{f(x_0+h, y_0+k) - f(x_0+h,y_0) - f(x_0,y_0 + k) + f(x_0,y_0)}{hk} - \frac{\partial^2 f}{\partial x \partial y}(x_0, y_0) \right| < \epsilon$$
Proof
Let $\epsilon_1, \epsilon_2, \epsilon_3, \epsilon > 0$, $0 < |h| < \delta$, and $0 < |k| < \delta$.\
Let $\varphi(x) = f(x, y_0+k) - f(x,y_0)$, where $x\in \mathbb{R}$.\
$\varphi$ is differentiable since 2nd-order partials are continuous, so we can apply the Mean Value Theorem for $x\in [x_0, x_0+h]$ to select a point $x_1$ in the open interval $(x_0, x_0+h)$ such that $$\frac{\varphi(x_0 + h) - \varphi(x_0)}{h} = \varphi'(x_1) = \frac{\partial f}{\partial x}(x_1, y_0 + k) - \frac{\partial f}{\partial x}(x_1, y_0)$$
By continuity, since $|h|<\delta$, $$\left| \frac{\varphi(x_0 + h) - \varphi(x_0)}{h} - \varphi'(x_1) \right| < \epsilon_1$$
We can define another auxiliary function $\alpha(y) = \frac{\partial f}{\partial x}(x_1, y)$, $y\in \mathbb{R}$\
Again, we can apply the mean value theorem to $y\in [y_0, y_0 + k]$ to find a point $y_1$ in that open interval such that $$\frac{\alpha(y_0 + k) - \alpha(y_0)}{k} = \frac{\partial^2 f}{\partial y \partial x}(x_1, y_1)$$
By continuity, since $|k|<\delta$, $$\left| \frac{\alpha(y_0 + k) - \alpha(y_0)}{k} - \frac{\partial^2 f}{\partial y \partial x}(x_1, y_1) \right| < \epsilon_2$$
Also by continuity, of 2nd-order derivatives, $$\left|\frac{\partial^2 f}{\partial y \partial x}(x_1, y_1) - \frac{\partial^2 f}{\partial y \partial x}(x_1, y_1)\right| < \epsilon_3$$
Note that $$\frac{\alpha(y_0 + k) - \alpha(y_0)}{k} = \frac{f(x_0+h, y_0+k) - f(x_0+h,y_0) - f(x_0,y_0 + k) + f(x_0,y_0)}{hk}$$
If $\epsilon = \max\{\epsilon_1, \epsilon_2, \epsilon_3\}$, then the desired result is proven by the last equation.