I have to show that following discretization of a linear differential operator satisfies the equation $\Delta_h u(x) = \Delta u(x) + \frac{h^2}{12} \Delta^2 u(x) + O(h^4)$
$$\Delta_h u = \frac{1}{6h^2} \begin{bmatrix} 1 & 4 & 1 \\ 4 & -20 & 4 \\ 1 & 4 & 1 \end{bmatrix}_* u$$
Now the conventional Laplace operator $\Delta$ results in a $3\times 3$ kernel, but the operator $\Delta^2$ has fourth order derivatives and therefore has a kernel size of at least $5\times 5$ is that correct? (My idea was $\Delta^2 u = \Delta (u_{xx} + u_{yy}) = u_{xxxx}+2u_{xxyy}+u_{yyyy}$, is that correct?)
This could be discretisized as follows (step size will always be $h$ here):
$$u_{xx} = \frac{1}{h^2}(u(x+h)-2u(x)+u(x-h)+O(h^2)$$
$$u_{xxxx} = \frac{1}{h^2}(u_{xx}(x+h)-2u_{xx}(x)+u_{xx}(x-h)+O(h^2)$$
Therefore
$$u_{xxxx} = \frac{1}{h^4}[1 ,-4 , 6 ,-4 ,1]_* u + O(h^2)$$
analogous for $u_{yyyy}$
Then I calculated the discretization for $u_{xxyy}$ which resulted in
$$u_{xxyy} = \frac{1}{h^2}\begin{bmatrix} 1 & -2 & 1 \\ -2 & 4 & -2 \\ 1 & -2 & 1 \end{bmatrix}_* u + O(h^2)$$
So
$$\Delta^2 =u_{xxxx}+2u_{xxyy}+u_{yyyy} = \frac{1}{h^2}\begin{bmatrix} &&1&& \\&2 & -8 & 2 &\\1& -8 & 20 & 8 & 1 \\ & 2 & -8 & 2 & \\ && 1 &&\end{bmatrix}_*+O(h^2)$$
But now I have a $5\times 5$ kernel and I want the final kernel to be only $3\times 3$. So my idea was to eliminating the 1 outside by dividing the laplace operator up into a discretization with stepsize $2h$ and one with $h$.
I found that
$$\Delta u = \Delta_{2h} u + O(h^2) = \frac{1}{4h^2}\begin{bmatrix} &&1&& \\& & 0 & &\\1& 0 & -4 & 0 & 1 \\ && 0 && \\ && 1 &&\end{bmatrix}_* u +O(h^2) $$
and
$$\Delta u = \Delta_{h} u + O(h^2) = \frac{1}{h^2}\begin{bmatrix} &1& \\1 & -4 & 1 \\& 1 &\end{bmatrix}_* u +O(h^2) $$
So I could write it as
$$\Delta_h u= \frac 8 6 \Delta_h u - \frac 2 6 \Delta_{2h} u + \frac{h^2}{12} \Delta_h^2 u$$
But now I get here an error of $O(h^4)$ for the expression $\frac{h^2}{12} \Delta_h^2 u$ but only $O(h^2)$ for the expressions $\frac 8 6 \Delta_h u - \frac 2 6 \Delta_{2h}$ is that correct or did I do something wrong? If everything was correct and my result is really only $O(h^2)$ how could I achieve $O(h^4)$?
As BaronVT commented, $$ g''(t) = \frac{1}{h^2}(g(t+h) - 2g(t) + g(t-h)) + \frac{h^2}{12}g^{(4)}(t) + O(h^4). $$ Now apply this to $g(t)=u(x+t,y+t)$ at $t=0$ to get the first diagonal. Here $g''(0)=u_{xx}+2u_{xy}+u_{yy}$ and $g^{(4)}(0)=u_{xxxx}+4u_{xxxy}+6u_{xxyy}+4u_{xyyy}+u_{yyyy}$.
The same for the other diagonal $g(t)=u(x+t,y-t)$ has derivatives $g''(0)=u_{xx}-2u_{xy}+u_{yy}$ and $g^{(4)}(0)=u_{xxxx}-4u_{xxxy}+6u_{xxyy}-4u_{xyyy}+u_{yyyy}$ so that in the sum of both diagonals the "odd" terms cancel.
Now repeat the same for the axes $g(t)=u(x+t,y)$ and $g(t)=u(x,y+t)$ and sum everything up to hopefully get the claimed result.