= Solution
Let all spatial <derivatives> be evaluated at $x=mh$, $t=nk$. Expand the two adjacent flux contributions symmetrically. The bracket in the scheme is
$$
S=a(x-h/2)[u(x-h)-u(x)]+a(x+h/2)[u(x+h)-u(x)]
=h^2(au_x)_x+O(h^4).
$$
For clarity, the $h^2$ terms are $au_{xx}+a'u_x$; odd orders cancel under $h\mapsto-h$. Under four spatial <derivatives> the next coefficient is $au_{xxxx}/12+a'u_{xxx}/6+a''u_{xx}/8+a^{(3)}u_x/24$. Time expansion gives $u(x,t+k)=u+k u_t+O(k^2)$. Since $u_t=(au_x)_x$ and $k=\mu h^2$,
$$
u(x,t+k)-u(x,t)-\mu S=O(k^2)+O(\mu h^4)=\boxed{O(h^4)}.
$$
This stronger smooth-data residual implies the requested \b[$O(h^3)$] consistency bound. Even a less detailed Taylor expansion with an $O(h^3)$ spatial remainder would suffice for the printed order. Uniform bounds on the needed <derivatives> over the region under consideration make the estimate uniform in the interior mesh index; the zero boundary values enter the scheme exactly.
Back to article page
