Interesting question. I too would like to know if this is mentioned anywhere explicitly. My guess is maybe this is a little too hard and not too useful to be used/taught. Perhaps someone can give a proper reference.
My calculations (see code below) suggest, if I am not wrong, that for the equation
$$ u_t = \mathcal{L}u, \qquad \mathcal{L}u = \nabla\cdot(D\nabla u), $$
there is such a nine-point stencil $L_9$. I am going to assume the central point of the stencil is $(0,0)$ and examine that point.
It's a quick and dirty calculation, which is why it (unusually) uses derivatives of $D$ in the stencil coefficients.
Here are its coefficients (sorry about the formatting; also the matrices were rotated $90^\circ$ clockwise by accident (sorry!), but they are a pain to retype):
$$
\begin{aligned}
&
\frac{D}{6h^2}
\begin{pmatrix}
1 & 4 & 1 \\
4 & -20 & 4 \\
1 & 4 & 1
\end{pmatrix}
\\&+
\frac{1}{12h}
\begin{pmatrix}
-D_y-D_x & -4D_x & D_y-D_x \\
-4D_y & 0 & 4D_y \\
-D_y + D_x & 4D_x & D_y+D_x
\end{pmatrix}
\\&+
\frac{1}{24D}
\begin{pmatrix}
-D_x D_y & -4 D_x^2-2 D_y^2 & D_x D_y \\
-2 D_x^2-4 D_y^2 & 12 D_x^2+12 D_y^2 & -2 D_x^2-4 D_y^2 \\
D_x D_y & -4 D_x^2-2 D_y^2 & -D_x D_y
\end{pmatrix}
\\&+
\frac{1}{12}
\begin{pmatrix}
D_{{xy}} & 3 D_{{xx}}+D_{{yy}} & -D_{{xy}} \\
D_{{xx}}+3 D_{{yy}} & -8 D_{{xx}}-8 D_{{yy}} & D_{{xx}}+3 D_{{yy}} \\
-D_{{xy}} & 3 D_{{xx}}+D_{{yy}} & D_{{xy}}
\end{pmatrix}
\\&+
\frac{h}{24D}
\begin{pmatrix}
0 & D_x D_{{xx}}+D_{{xy}} D_y & 0 \\
D_x D_{{xy}}+D_y D_{{yy}} & 0 & -D_x D_{{xy}}-D_y D_{{yy}} \\
0 & -D_x D_{{xx}}-D_{{xy}} D_y & 0
\end{pmatrix}
\\&+
\frac{h}{24}
\begin{pmatrix}
0 & -D_{xyy} - D_{xxx} & 0 \\
-D_{yyy} - D_{xxy} & 0 & D_{yyy} + D_{xxy} \\
0 & D_{xyy} + D_{xxx} & 0
\end{pmatrix}
\end{aligned}
$$
It's much like the original stencil, but with some corrections. Unfortunately, they look a bit awkward to me, so I can't say much about them beyond saying that they satisfy the equations I need them to satisfy.
I calculated the coefficient by brute force: by writing down the linear equations satisfied by the stencil coefficients, such that the approximation at the point $(0,0)$ would satisfy
$$ L_9u = \mathcal{L}u + \beta h^2\nabla^2\mathcal{L}u + \gamma h^2 \mathcal{L}\mathcal{L}u + O(h^4). $$
The point is that this kind of method of deferred corrections only works if in the equation $\mathcal{L}u = \cdots$ the first error term of order $h^2$ looks like $h^2\mathcal{M}\mathcal{L}u$ for some suitable operator $\mathcal{M} = \beta \nabla^2 + \gamma \mathcal{L}$.
I then calculated the error term coefficients
$$ \beta=\frac{1}{6}, \qquad \gamma = \frac{-1}{12D(0,0)}. $$
The corresponding modified equation after space discretization is
$$ \left(1 + \frac{h^2}6\nabla^2 - \frac{h^2}{12D(0,0)}\mathcal{L} \right) u_t = L_9 u. $$
Here is the code that I used:
(* VarCoefNinePointStencil.nb *)
Clear[ncdiff, ncdiff2]
ncdiff[d_, u_] :=
D[d[x, y] D[u[x, y], x], x] + D[d[x, y] D[u[x, y], y], y]
ncdiff2[d_, u_] := D[d[x, y] D[u, x], x] + D[d[x, y] D[u, y], y]
Module[{L9, Lapprox, trunc, conds, sol, target, errTerm1, errTerm2,
uvars, remainder, remCond, remSol},
L9 = Table[Subscript[\[Alpha], i, j], {i, -1, 1}, {j, -1, 1}];
Lapprox = Normal@Series[
Sum[L9[[i + 2, j + 2]] u[i h, j h], {i, -1, 1}, {j, -1, 1}]
, {h, 0, 4}];
uvars = {Derivative[a_, b_][u][0, 0] :> Subscript[u, a, b],
u[0, 0] -> Subscript[u, 0, 0]};
Lapprox = Lapprox /. uvars;
errTerm1 = ncdiff[d, u];
errTerm2 = ncdiff[d, u];
errTerm1 = D[errTerm1, x, x] + D[errTerm1, y, y];
errTerm2 = ncdiff2[d, errTerm2];
target =
ncdiff[d, u] + h^2 \[Beta] errTerm1 +
h^2 \[Gamma] errTerm2 /. {x -> 0, y -> 0};
target = target /. uvars;
trunc = Lapprox - target /. {x -> 0, y -> 0};
trunc = trunc /. uvars;
conds = Flatten@
Table[Coefficient[trunc, Subscript[u, a, b]], {a, 0, 2}, {b, 0, 2}];
Print[conds // Column];
sol = Solve[Thread[conds == 0], Flatten@L9];
If[Head[sol] === Solve || sol === {}, Print["No solution."]; Abort[]];
sol = sol[[1]];
Print[sol];
Print[Lapprox - target /. sol // Expand // Collect[#, h] &];
remainder = Coefficient[Lapprox - target /. sol, h^2];
remCond =
Coefficient[remainder, #] & /@ {Subscript[u, 4, 0], Subscript[u, 3,
0]};
remSol = Solve[Thread[0 == remCond], {\[Beta], \[Gamma]}];
Print[remSol];
Print[(1/(24 h^2)) MatrixForm[
24 h^2 L9 //. Join[sol, remSol[[1]]] // Expand]];
Print[Lapprox - target //. Join[sol, remSol[[1]]] //
Collect[#, h, Factor] &];
Collect[remainder,
Union@Cases[remainder, Subscript[u, _, _], \[Infinity]], Factor]
]
