There are two main ways to see finite differences. One way is to see it as taking linear combinations of the Taylor expansions of the points and choosing the coefficients so that you get a desired order of consistency. The other way (which is more or less equivalent) is to take the derivative of the interpolating polynomial of the points. Thus having your nine points you can compute the Taylor expansions:
\begin{align}
u(x+a,y+b) &= u(x,y) + (a\partial_x + b\partial_y) u(x,y) + \frac{1}{2!} (a\partial_x + b\partial_y)^2 u(x,y) + O(\max(a,b)^3)\\
& = \sum_{k=1}^{n} \frac{1}{k!} (a\partial_x + b\partial_y)^k u(x,y) + O(\max(a,b)^{n+1})
\end{align}
Now you introduce one coefficient per point and you take the linear combinations, in your case:
$$S = \sum_{i=-1}^{j=1}\sum_{i=-1}^{j=1} \alpha_{i,j} u(x+a,y+b) = \frac{\partial^2 u}{\partial x \partial y} + O(h^p).$$
Expanding this you get:
\begin{align}
S &= \left(\sum_{i=-1}^{1}\sum_{j=-1}^1 \alpha_{i,j}\right)u(x,y)
\\
&+\left(\sum_{i=-1}^{1}\sum_{j=-1}^1 (ih)\alpha_{i,j}\right)\partial_{x} u(x,y) + \left(\sum_{i=-1}^{1}\sum_{j=-1}^1 (jh)\alpha_{i,j}\right)\partial_{y} u(x,y)
\\
&+\left(\sum_{i=-1}^1\sum_{j=-1}^1 \frac{(ih)^2}{2} \alpha_{i,j}\right)\partial^2_x u(x,y)
+ \left(\sum_{i=-1}^1\sum_{j=-1}^1 2\frac{(ih)(jh)}{2} \alpha_{i,j}\right)\partial_x\partial_y u(x,y)
\\
& + \left(\sum_{i=-1}^1\sum_{j=-1}^1 \frac{(jh)^2}{2} \alpha_{i,j}\right)\partial^2_y u(x,y) +
\left(\sum_{i=-1}^1\sum_{j=-1}^1 \frac{(ih)^3}{3!} \alpha_{i,j}\right)\partial^3_x u(x,y) \\
&+
\left(\sum_{i=-1}^1\sum_{j=-1}^1 3\frac{(ih)^2(jh)}{3!} \alpha_{i,j}\right)\partial^2_x\partial_y u(x,y)
+\left(\sum_{i=-1}^1\sum_{j=-1}^1 3\frac{(ih)(jh)^2}{3!} \alpha_{i,j}\right)\partial_x\partial^2_y u(x,y) \\
&+\left(\sum_{i=-1}^1\sum_{j=-1}^1 \frac{(jh)^3}{3!} \alpha_{i,j}\right)\partial^3_y u(x,y) + O(h^4).
\end{align}
Since you want only $\partial_x \partial_y u$ all coefficients in front of the other terms should be zero:
$$\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\
-h & 0 & h & -h & 0 & h & -h & 0 & h \\
-h & -h & -h & 0 & 0 & 0 & h & h & h \\
\frac{h^2}{2} & 0 & \frac{h^2}{2} & \frac{h^2}{2} & 0 & \frac{h^2}{2} & \frac{h^2}{2} & 0 & \frac{h^2}{2}\\
2\frac{h^2}{2} & 0 & 2\frac{h^2}{2} & 0 & 0 & 0 & 2\frac{h^2}{2} & 0 & 2\frac{h^2}{2} \\
\frac{h^2}{2} & \frac{h^2}{2} & \frac{h^2}{2} & 0 & 0 & 0 & \frac{h^2}{2} & \frac{h^2}{2} & \frac{h^2}{2} \\
-\frac{h^3}{6} & 0 & \frac{h^3}{6} & -\frac{h^3}{6} & 0 & \frac{h^3}{6} & -\frac{h^3}{6} & 0 & \frac{h^3}{6} \\
-3\frac{h^3}{6} & 0 & -3\frac{h^3}{6} & 0 & 0 & 0 & 3\frac{h^3}{6} & 0 & 3\frac{h^3}{6} \\
-3\frac{h^3}{6} & 0 & 3\frac{h^3}{6} & 0 & 0 & 0 & -3\frac{h^3}{6} & 0 & 3\frac{h^3}{6} \\
-\frac{h^3}{6} & -\frac{h^3}{6} & -\frac{h^3}{6} & 0 & 0 & 0 & \frac{h^3}{6} & \frac{h^3}{6} & \frac{h^3}{6}
\end{bmatrix}
\begin{bmatrix} \alpha_{-1,-1} \\ \alpha_{0,-1} \\ \alpha_{1,-1} \\ \alpha_{-1,0} \\ \alpha_{0,0} \\ \alpha_{1,0} \\ \alpha_{-1,1} \\ \alpha_{0,1} \\ \alpha_{1,1} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0\end{bmatrix} \begin{matrix} (u) \\ (\partial_x u) \\ (\partial_y u) \\ (\partial^x_2 u) \\ (\partial_x\partial_y) \\ (\partial_y^2) \\ (\partial_x^3) \\ (\partial_x^2\partial_y) \\ (\partial_x \partial^2_y) \\ (\partial_y^3)\end{matrix}.$$
You will notice that you have 10 equations and 9 unknowns, so this system is overdetermined. This is not surprising as the interpolation problem for a cubic bivariate polynomial requires 10 points. To get a solution you likely have to give up on one of the equations. This, however, is not very nice in terms of symmetry. What I would do is that I would solve the system by dropping the equations for $\partial^2_x\partial_y$ and $\partial_x \partial^2_y$ (provided you care about consistency along the main axes more), and then I would use the one free variable to find the best approximation that minimizes the coefficients in front of $\partial^2_x\partial_y$ and $\partial_x \partial^2_y$.If you care more about the diagonal than the main axes then I would suggest throwing away the equations for $\partial^3_x$ and $\partial^3_y$ and instead keeping $\partial^2_x\partial_y$ and $\partial_x \partial^2_y$. Then you can minimize the coefficients in front of $\partial^3_x$ and $\partial^3_y$ in order to get a unique solution.
