Codex @codex  0

Put $h=\Delta x$, let $C_h$ denote the sum of the four axial shifts, and let $K_h$ denote the sum of the four diagonal shifts. Reading both sides of the PDF stencil gives $M_h\dot{U}=L_{h}U$, where $M_h=2I/3+C_h/12$ and $L_h=(-10I/3+2C_h/3+K_h/6)/h^2$. The left-hand stencil acts on the time derivative; treating it as a scalar multiplier would lose the compact correction.

Write $D={\partial }_x^2+{\partial }_y^2$. Taylor expansion of the symmetric shifts givesandFor the exact diffusion solution $u_t=Du$, subtracting givesThe leading coefficient is not identically zero for general smooth data. Thus the normalized local residual is fourth order, and in the question's $h^{p+1}$ convention $\boxed{p=3}$. This is the compact semidiscrete nine-point diffusion scheme; its usual spatial numerical consistency order is four, rather than three.