The Lagrange interpolation on $\mathbb{R}^2$ of $f\in\mathcal{C}^{2}([a,b]\times[c,d])$ is defined as: $$f(x,y)\approx\sum_{m=0}^M \sum_{n=0}^{N}L_{n,m}(x,y) f(x_n,y_m),\qquad N,M\in\mathbb{N},$$ where $\{x_n\}$ is a partition of $[a,b]$, $\{y_m\}$ is a partition of $[c,d]$, and $$L_{n,m}(x,y)=\prod_{j=0,j\neq n}^N\dfrac{x-x_j}{x_{j+1}-x_j}\prod_{l=0,l\neq m}^M\dfrac{y-y_l}{y_{l+1}-y_l}.$$
If we consider the particular case of the Lagrange interpolation on $\mathbb{R}$, the error for $x\in[a,b]$ is $$E(x)=f(x)-H(x)=\dfrac{f''(\xi)}{2}\prod_{n=0}^N(x-x_n),\qquad \xi\in(a,b),$$ where $$H(x)=\sum_{n=0}^{M}f(x_n)\prod_{j=0,j\neq n}^N\dfrac{x-x_n}{x_{n+1}-x_n},\qquad N\in\mathbb{N}.$$
Then, I am wondering if there is a formula for the two dimensional case.