The ordinary difference quotient is the simplest formula for univariate numerical differentiation. Polynomial interpolation is one of the basic theories for those simple formulas such as forward difference $$f'(x)\approx \frac{f(x+h)-f(x)}{h},$$ backward difference $$f'(x)\approx \frac{f(x)-f(x-h)}{h},$$ and centered difference $$f'(x)\approx \frac{f(x+h)-f(x-h)}{2h}.$$ We use the derivative of an interpolation polynomial to approximate the derivative of $f(x)$, so that we can obtain these formulas. In the bivariate case, using the Lagrange interpolation polynomial on the grid-points, for example, the nine sample points are $$(x_0-h_1,y_0+h_2) , (x_0,y_0+h_2), (x_0+h_1,y_0+h_2)$$ $$(x_0-h_1,y_0) , (x_0,y_0), (x_0+h_1,y_0)$$ $$(x_0-h_1,y_0-h_2) , (x_0,y_0-h_2), (x_0+h_1,y_0-h_2)$$ we can obtain the numerical differentiation formulas \begin{eqnarray*} % \nonumber to remove numbering (before each equation) \frac{\partial}{\partial x}f(x_0,y_0) &=& \frac{f(x_0+h_1,y_0)-f(x_0-h_1,y_0)}{2h_1}, \\ \frac{\partial}{\partial y}f(x_0,y_0) &=& \frac{f(x_0,y_0+h_2)-f(x_0,y_0-h_2)}{2h_2}. \end{eqnarray*} It is very hard to obtain the numerical differentiation formulas when the points are not grid shape, such as the Delaunay triangulation. In this paper, based on multivariate polynomial interpolation, we give a numerical differentiation algorithm to compute the differential values of multivariate functions at a single point directly.