What is the remainder for a taylor series of two variables. I know that for a function of one variable $f(x)$, the Taylor expansion is $$f(x)=f(x_0)+f'(x_0)(x-x_0)+...+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x,x_0) $$
where 
\begin{align*}
R_n(x,x_0)&=\int^x_{x_0}\frac{(x-t)^n}{n!} f^{(n+1)}(t)dt \\
&=f^{(n+1)}(\xi)\int^x_{x_0} \frac{(x-t)^n}{n!}dt \\
&=f^{(n+1)}(\xi) \frac{(x-x_0)^{n+1}}{(n+1)!}.
\end{align*}
My question is for $f(x,y)$ with a Taylor expansion
$$f(x)=f(x_0,y_0)+f_x(x_0,y_0)(x-x_0)+f_y(x_0,y_0)(y-y_0)+...+R_n(x,x_0,y,y_0) $$ 
What is $R_n(x,x_0,y,y_0)$?
 A: We consider  $x,x_0,y,y_0\in\mathbb{R}$ and define intervals $I,J$ with
\begin{align*}
&I=[\mathrm{min}(x,x_0),\mathrm{max}(x,x_0)]\qquad\text{and}\qquad  J=[\mathrm{min}(y,y_0),\mathrm{max}(y,y_0)]
\end{align*}

Let  $f:I\times  J\rightarrow\mathbb{R}$ be an $(n+1)$-times continuously differentiable function in $x$  and $y$. There exist values   $\xi\in I^\circ,\eta\in  J^\circ$, so that following is valid
  \begin{align*}
f(x,y)&=f(x_0,y_0)+f_x(x_0,y_0)(x-x_0)+f_y(x_0,y_0)(y-y_0)\\
&\qquad+\frac{1}{2}\left(f_{xx}(x_0,y_0)(x-x_0)^2+2f_{xy}(x_0,y_0)(x-x_0)(y-y_0)\right.\\
&\qquad\qquad\quad\left.+f_{yy}(x_0,y_0)(y-y_0)^2\right)\\
&\qquad\,\,\,\vdots\\
&\qquad+\frac{1}{n!}\left(f_{x^n}(x_0,y_0)(x-x_0)^n+\binom{n}{1}f_{x^{n-1}y}(x_0,y_0)(x-x_0)^{n-1}(y-y_0)\right.\\
&\qquad\qquad\quad+\cdots+f_{y^n}(x_0,y_0)(y-y_0)^n\Big)\\
&\qquad+R_n(x,y)
\end{align*}
  The Lagrange remainder term $R_n(x,y)$  is given as
  \begin{align*}
\color{blue}{R_n(x,y)}&\color{blue}{=\frac{1}{(n+1)!}\left(f_{x^{n+1}}(\xi,\eta)(x-x_0)^{n+1}+\binom{n+1}{1}f_{x^{n-1}y}(\xi,\eta)(x-x_0)^n(y-y_0)\right.}\\
&\qquad\qquad\qquad\quad\color{blue}{+\cdots+f_{y^{n+1}}(\xi,\eta)(y-y_0)^{n+1}\Big)}
\end{align*}

A more compact notation can be given with the  total derivative. We obtain
\begin{align*}
df&=(x-x_0)f_x+(y-y_0)f_y,\\
d^2f&=((x-x_0)f_x+(y-y_0)f_y)^{(2)}\\
&=((x-x_0)^2f_{xx}+2 (x-x_0)(y-y_0)  f_{xy}+(y-y_0)^2f_{yy}\\
&\,\,\,\vdots\\
d^nf&=((x-x_0)f_x+(y-y_0)f_y)^{(n)}\\
&=(x-x_0)^nf_{x^n}+\binom{n}{1}(x-x_0)^{n-1}(y-y_0)f_{x^{n-1}y}\\
&\qquad+\cdots+(y-y_0)^nf_{y^n}
\end{align*}

and we can  write
  \begin{align*}
f(x,y)=f(x_0,y_0)+df(x_0,y_0)+\frac{1}{2}d^2f(x_0,y_0)+\cdots+\frac{1}{n!}d^nf(x_0,y_0)+R_n(x,y)
\end{align*}
  with
  \begin{align*}
\color{blue}{R_n(x,y)=\frac{1}{(n+1)!}d^{n+1}f(\xi,\eta)}
\end{align*}

Note:  A  derivation  similar to this answer can be found in section 6.3  of  Introduction to Calculus and Analysis II by Richard Courant.
