$\mathbf {Theorem}$: if $\ f$ is a function of $\ x$ and $\ y$, where $\ f_{x}$ and $\ f_{y}$ are continuous in an open region $\ R$, then $\ f$ is differentiable on $\ R$.
Proof:Let $\ S$ be a surface defined by $\ z=f(x,y)$, where $\ f$, $\ f_{x}$, and $\ f_{y}$ are continuous at $\ (x,y)$.
Let $\Delta z$ be the total change in $\ f$ from two points on the surface $\ S$ , let $\Delta z_{x}$ be the change in $\ f$ when $\ y$ is held constant and let $\Delta z_{y}$ be the change in $\ f$ when x is held constant. $$\Delta z=f(x+\Delta x, y+\Delta y)-f(x,y)$$ $$\Delta z= [f(x+\Delta x,y)-f(x,y)]+[f(x+\Delta x, y+\Delta y)-f(x+\Delta x,y)]$$ $$\Delta z= \Delta z_{x}+\Delta z_{y}$$
By the Mean Value Theorem there is some $\ x_{1}$ between $\ x$ and $\ x+\Delta x$ and some $\ y_{1}$ between $\ y$ and $\ y+\Delta y$ such that $$\Delta z_{x}=f_{x}(x_{1},y)\Delta x$$ $$\Delta z_{y}=f_{y}(x+\Delta x,y_{1})\Delta y$$
If you define $\epsilon_{1}$ and $\epsilon_{2}$ as $$\epsilon_{1}=f_{x}(x_{1},y)-f_{x}(x,y)$$ $$\epsilon_{2}=f_{y}(x+\Delta x, y_{1})-f_{y}(x,y)$$
It follows that $$\Delta z=\Delta z_{x}+\Delta z_{y}=[f_{x}(x,y)\Delta x + f_{y}(x,y)\Delta y]+\epsilon_{1}\Delta x + \epsilon_{2}\Delta y$$
By the continuity of $\ f_{x}$ and $\ f_{y}$ and the fact that $\ x\le x_{1}\le x+\Delta x$, and $\ y\le y_{1}\le y+\Delta y$, it follows that $$\epsilon_{1}\rightarrow 0\ \ \text{and}\ \ \epsilon_{2}\rightarrow 0\ \ \text{as}$$ $$\Delta x\rightarrow 0\ \ \text{and}\ \ \Delta y\rightarrow 0$$
Therefore, by definition, $\ f$ is differentiable.
My question is, why do the continuity of the partial derivatives $\ f_{x}$ and $\ f_{y}$ and the intervals that $\ x_{1}$ and $\ y_{1}$ lie on imply that $\ f$ is differentiable?
By definition, $f$ is differentiable at $(x,y)$ with derivative $(f_x(\vec x), f_y(\vec x))$ iff $$ \lim_{(\Delta x, \Delta y) \to (0,0)} \frac{|f(x +\Delta x, y + \Delta y) - f( x, y) - f_x(x, y)\Delta x - f_y(x , y) \Delta y|}{|(\Delta x, \Delta y)|} = 0$$
But by your calculation, $$ \frac{|f(x +\Delta x, y + \Delta y) - f( x, y) - f_x(x, y)\Delta x - f_y(x , y) \Delta y|}{|(\Delta x, \Delta y)|} = \frac{|\epsilon_1 \Delta x + \epsilon_2 \Delta y|}{\sqrt{(\Delta x)^2 + (\Delta y)^2}}\leq |\epsilon_1 | + |\epsilon_2|,$$ where $\epsilon_1$ and $\epsilon_2$ are as defined in your original post.
So it suffices to show that $ \lim_{(\Delta x, \Delta y) \to (0,0)} \epsilon_1 = \lim_{(\Delta x, \Delta y) \to (0,0)} \epsilon_2 = 0.$
The equality $$ \lim_{(\Delta x, \Delta y) \to (0,0)} \epsilon_1 = \lim_{(\Delta x, \Delta y) \to (0,0)} (f_x(x_1, y) - f_x (x, y)) = 0$$ holds by composition of limits. First, since $x \leq x_1 \leq x + \Delta x$, we have $\lim_{\Delta x \to 0} x_1 = x$. Second, $\lim_{x_1 \to x} f_x(x_1, y) = f_x(x, y)$ holds by the continuity of $f_x$ at $(x,y)$. Composing the two limits gives the desired result.
The equality $$ \lim_{(\Delta x, \Delta y) \to (0,0)} \epsilon_2 = \lim_{(\Delta x, \Delta y) \to (0,0)} (f_x(x + \Delta x, y_1) - f_x (x, y)) = 0$$ follows in the same way. Let $x_2 = x + \Delta x$. First, observe that $\lim_{\Delta x \to 0}x_2 = x$, and observe that $\lim_{\Delta y \to 0}y_1 = y$ since $y \leq y_1 \leq y + \Delta y$. Second, observe that $\lim_{(x_2, y_1) \to (x, y)} f_y(x_2,y_1) = f_y(x, y)$, by the continuity of $f_y$ at $(x,y)$. Then compose the two limits to get the desired result.