In 2D you would need a bi-quadratic interpolation function such as $(a_0 + a_1 x + a_2 x^2)(b_0 + b_1 y + b_2 y^2)$ giving rise to 9 undetermined coefficients. Thus you need 9 points where the values are known then you can determine these coefficients.
At the least you need incomplete polynomial up to 2nd order for a 2nd order estimate, i.e., $\{1,x,y,x^2,y^2,xy\}$ with 6 coefficients, thus you need 6 points.
With 5 points and 6 coefficients you will have an under determined linear least squares problem. You can convert it to an over determined system by adding more points, especially if you are on, say, a finite difference grid. You can find details in any numerical analysis text book. See also this. I give some basics below.
Your goal is to find the polynomial function coefficients - the function itself is linear in the coefficients. Now express the sum of the square of error between interpolated values $v^I$ and actual values $v^A$ at $n>6$ points as:
\begin{align*}
E = \sum_{i=1}^n(v^I(a_p, \mathbf x_i) - v^A_i)^2 \qquad& v^I \text{ is linear in } a_p \text{ evaluated at } \mathbf x_i \\
= \sum_{i=1}^n(v^I_p(\mathbf x_i) a_p - v^A_i)^2\qquad& \text{p}^{th} (1...6) \text{ coefficient of } v^I \text{ at } \mathbf x_i \\
{{\partial E}\over{\partial a_k}} = \sum_{i=1}^n 2(v^I_p(\mathbf x_i) a_p - v^A_i) v^I_p(\mathbf x_i)\delta_{pk} = 0\qquad& \text{minimize w.r.t. } a_k \\
\implies \sum_{i=1}^n v^I_k(\mathbf x_i) v^I_p(\mathbf x_i) a_p = \sum_{i=1}^n v^I_k(\mathbf x_i) v^A_i\qquad&
\end{align*}
Clearly if you have an over determined system, you can also go with the approach of bi-quadratic functions with 9 coefficients (on a 3x3 grid) at each point.
