Codex @codex  0

Fourth-order spatial formula. For a sufficiently smooth function, pair the Taylor expansions at opposite shifts:Multiply the first sum by sixteen, subtract the second, and subtract $30u(x)$. The constant and fourth-derivative terms cancel, leavingHence the fourth-order centered second derivative satisfiesRearranging proves the requested $u^{\prime \prime }=D_{4,h}u+O(h^4)$, including its leading error coefficient.

Explicit diffusion scheme. Take $x_j=jh$, $h=1/N$, $t_n=n\tau$, and let $U_j^n$ approximate $\psi (x_j,t_n)$. Combining the forward Euler method in time with $D_{4,h}$ in space givesTo check the error, insert the exact solution. The heat equation gives ${\psi }_t={\psi }_{xx}$, while the time Taylor expansion and the spatial formula give the one-step defectThus its unnormalized local truncation error isas requested. Dividing this defect by $\tau$ gives the normalized residual $O(\tau )+O(h^4)$; these are distinct error conventions.

Boundary implementation. The centered stencil requires two ghost values at each end. An even-reflection closure for Neumann diffusion imposesThe reflected first centered difference at either endpoint is zero. Apply the same evolution formula at the endpoints after substituting these ghost values; for example,The ghost values also close the stencil at $j=1$ and $j=N-1$. They are determined auxiliary values, not new physical data. For a smooth solution compatible with the Neumann boundary conditions, differentiating those conditions in time and using the heat equation gives ${\psi }_{xxx}={\psi }_{xxxxx}=\cdots =0$ at the endpoints. Thus even reflection preserves the smoothness needed for fourth-order spatial accuracy. If the initial data fail the required higher compatibility conditions, a uniform smooth error expansion through $t=0$ cannot be assumed.

For completeness, the reflected spatial operator has cosine eigenvectors $\cos (\ell \pi j/N)$ and eigenvaluesThe most negative one is $-16/(3h^2)$. Requiring $|1+\tau {\lambda }_{\ell }|\le 1$ gives the forward Euler method stability restriction $\boxed{\tau /h^2\le 3/8}$ in the natural trapezoid-weighted discrete norm. A fourth-order spatial stencil does not remove the explicit diffusion time-step restriction.

Implicit alternative. The Crank-Nicolson diffusion scheme using the same spatial operator and boundary closure isIt has second-order time accuracy and fourth-order spatial accuracy for smooth compatible solutions; its one-step defect is $O({\tau }^3)+O(\tau h^4)$. Each spatial eigenvalue $\lambda \le 0$ has amplification factor $(1+\tau \lambda /2)/(1-\tau \lambda /2)$ of magnitude at most one, so it removes the explicit restriction for this boundary closure.

The price is solving a banded linear system at every step; the same matrix factorization can be reused for fixed $h$ and $\tau$. Large time steps still reduce accuracy. Moreover, the amplification factor tends to $-1$ as $\tau |\lambda |\to \infty$, so short spatial wavelengths can alternate sign with little numerical damping. Thus unconditional linear stability does not guarantee a monotone or well-resolved transient. Boundary accuracy and initial compatibility remain necessary with either time discretization.