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)]'$ involves the mid point values of $a(x)$ and has seven stencils. Is there a formula that involves even less computations / stencils?