fourth-order finite difference for $(a(x)u'(x))'$ Previously I asked here about constructing a symmetric matrix for doing finite difference for $(a(x)u'(x))'$ where the (diffusion) coefficient $a(x)$ is spatially varying. The answer provided there works for getting a second order accurate method. What about getting a fourth-order accurate method? If I follow the same idea and apply the following fourth-order accurate formula for first derivative
$ u'_i = \dfrac{u_{i-1} - 8u_{i-1/2} + 8u_{i+1/2} - u_{i+1}}{6\Delta x}$   (1)
in succession, then I end up getting a formula for $(a u')_i'$ which involves $u_{i-2}, u_{i-3/2}, u_{i-1}, u_{i-1/2}, u_i, u_{i+1/2}, u_{i+1}, u_{i+3/2}, u_2$. It involves the mid point values $u_{i+n/2}$ and it depends on nine neighbouring values of $u_i$, which doesn't sound right for fourth order accurate scheme. For the special case of $a(x)=1$ it also doesn't reduce to the formula
$ u'' = \dfrac{-u_{i-2} + 16u_{i-1} - 30u_i + 16u_{i+1} - u_{i+2}}{12\Delta x^2}$.   (2)
So what's wrong with applying (1) in succession? What's the correct approach to get a finite difference formula for $(au')'$ with fourth order accuracy?
Edit 1: After reading this post I managed to derive a symmetric four-order centered difference scheme for $[a(x)u'(x)]'$ by applying
$ u'(x) = \dfrac{u_{i-3/2} - 27u_{i-1/2} + 27u_{i+1/2} - u_{i+3/2}}{24\Delta x} + \mathcal{O}(\Delta x^4)$
twice in succession. The resulting formula for $[a(x)u'(x)]'$ is
$ (au')'_i = \dfrac{1}{576\Delta x^2} \left( a_{i-\frac{3}{2}}u_{i-3}
- 27(a_{i-\frac{3}{2}} + a_{i-\frac{1}{2}})u_{i-2}
+ (27a_{i-\frac{3}{2}} + 729a_{i-\frac{1}{2}}+27a_{i+\frac{1}{2}})u_{i-1}
- (a_{i-\frac{3}{2}} + 729a_{i-\frac{1}{2}} + 729a_{i+\frac{1}{2}} + a_{i+\frac{3}{2}})u_i
+ (27a_{i-\frac{1}{2}} + 729a_{i+\frac{1}{2}}+27a_{i+\frac{3}{2}})u_{i+1}
- 27(a_{i+\frac{1}{2}} + a_{i+\frac{3}{2}})u_{i+2}
+ a_{i+\frac{3}{2}}u_{i+3} \right)$
which involves the mid point values of $a(x)$ and has seven stencils. Is there a formula that involves even less computations / stencils(e.g. five stencils)?
 A: Let's use
$$
f'(x) = \frac{f_{i-3/2} - 27 f_{i-1/2} + 27 f_{i+1/2} - f_{i+3/2}}{24 \Delta x} + O(\Delta x^4)
$$
as an approximation for the outer derivative. To obtain a fourth-order approximation for the $(au_x)_x$ term we need to obtain $a u_x$ at each of the $x_{i-3/2}, x_{i-1/2}, x_{i+1/2}, x_{i+3/2}$ points with fourth-order accuracy.
Assuming that we know $a$ everywhere, the problem reduces to computing $u_x$ at those points. Taking our stencil to be $(x_{i-2}, x_{i-1}, x_{i}, x_{i+1}, x_{i+2})$ we can compute fourth degree interpolating polynomial $P(x)$ for $u(x)$:
$$u(x) = P(x) + O(\Delta x^5).$$
Differentiating that polynomial gives fourth-order approximation for $u_x$:
$$
u_x(x) = P'(x) + O(\Delta x^4).
$$
Omitting straightforward computations (which I've done using Mathics) we get the following approximations for $u_x(x_{i-3/2}), \dots, u_x(x_{i+3/2})$
$$
\begin{pmatrix}
u_x(x_{i-3/2})\\
u_x(x_{i-1/2})\\
u_x(x_{i+1/2})\\
u_x(x_{i+3/2})
\end{pmatrix}
= 
\frac{1}{24 \Delta x}
\begin{pmatrix}
 -22 & 17 & 9 & -5 & 1\\ 
1 & -27 & 27 & -1 & 0 \\
0 & 1 & -27 & 27 & -1\\
 -1 & 5 & -9 & -17 & 22
\end{pmatrix}
\begin{pmatrix}
u_{i-2}\\
u_{i+1}\\
u_{i}\\
u_{i+1}\\
u_{i+2}
\end{pmatrix}
+ O(\Delta x^4).
$$
Combining everything to a single formula we get
$$
(a u_x)_x = \frac{1}{24 \Delta x} \times \left[\\
a_{i-3/2} \frac{u_{i+2} - 5 u_{i+1} + 9 u_i + 17 u_{i-1} - 22 u_{i-2}}{24 \Delta x}
-27a_{i-1/2} \frac{-u_{i+1} + 27 u_i - 27 u_{i-1} + u_{i-2}}{24 \Delta x} +\\
+27a_{i+1/2} \frac{-u_{i+2} + 27 u_{i+1} - 27 u_{i} + u_{i-1}}{24 \Delta x} -a_{i+3/2} \frac{22u_{i+2} - 17 u_{i+1} - 9 u_i + 5 u_{i-1} - u_{i-2}}{24 \Delta x}\\
\right] + O(\Delta x^4).
$$
A: A bit of work with five-point Lagrange interpolants based on a uniform grid for $a$ and $u$ shows, with
\begin{align}
w_{-2} &:= \tfrac{1}{144}a_{-2} - \tfrac{1}{18} a_{-1} - \tfrac{1}{12}a_0 + \tfrac{1}{18} a_1 - \tfrac{1}{144} a_2, \\
w_{-1} &:= -\tfrac{1}{18}a_{-2} + \tfrac{4}{9} a_{-1} + \tfrac{4}{3} a_0 - \tfrac{4}{9}a_1 + \tfrac{1}{18}a_2 \\
w_0 &:= -\tfrac{5}{2}a_0, \\
w_{1} &:= \tfrac{1}{18}a_{-2} - \tfrac{4}{9}a_{-1} + \tfrac{4}{3}a_0 + \tfrac{4}{9}a_1 - \tfrac{1}{18}a_2, \\
w_{2} &:= -\tfrac{1}{144}a_{-2} + \tfrac{1}{18}a_{-1} - \tfrac{1}{12}a_0 - \tfrac{1}{18}a_1 + \tfrac{1}{144}a_2
\end{align}
that $(au')'(0) = (w_{-2}u_{-2}+ w_{-1}u_{-1} + w_0u_0 + w_1 u_1 + w_2 u_2)/h^2 + O(h^4)$. It typically won't be a symmetric rule.
ADDED EDIT
Since I was thinking in terms of using grid values I started with the Lagrange interpolant for $a$
\begin{align*}
        A(x) &= \\
        & \quad a_{-2}\,
\frac{x - x_{-1}}{x_{-2} - x_{-1}}
\frac{x - x_{0}}{x_{-2} - x_{0}}
\frac{x - x_{1}}{x_{-2} - x_{1}}
\frac{x - x_{2}}{x_{-2} - x_{2}} \\
        &+  a_{-1}\,
\frac{x - x_{-2}}{x_{-1} - x_{-2}}
\frac{x - x_{0}}{x_{-1} - x_{0}}
\frac{x - x_{1}}{x_{-1} - x_{1}}
\frac{x - x_{2}}{x_{-1} - x_{2}} \\
        &+ a_{0}\,
\frac{x - x_{-2}}{x_{0} - x_{-2}}
\frac{x - x_{-1}}{x_{0} - x_{-1}}
\frac{x - x_{1}}{x_{0} - x_{1}}
\frac{x - x_{2}}{x_{0} - x_{2}} \\
        &+ a_{1}\,
\frac{x - x_{-2}}{x_{1} - x_{-2}}
\frac{x - x_{-1}}{x_{1} - x_{-1}}
\frac{x - x_{0}}{x_{1} - x_{0}}
\frac{x - x_{2}}{x_{1} - x_{2}} \\
        &+ a_{2}\,
\frac{x - x_{-2}}{x_{2} - x_{-2}}
\frac{x - x_{-1}}{x_{2} - x_{-1}}
\frac{x - x_{0}}{x_{2} - x_{0}}
\frac{x - x_{1}}{x_{2} - x_{1}}
\end{align*}
and the corresponding $U$ for $u$. Setting $x_{k} = kh$ for a uniform mesh simplifies the expressions and makes the template symmetric about $x=0$. Computing $(AU')'$ and evaluating it at $x=0$ gives the result which can be verified to be $O(h^4)$. The symmetric template is key.
