Solution (source code)

= Solution

Use a smooth exact solution such as $u(x,t)=e^{-Dq^2t}\sin(qx)$ on a domain with matching boundary conditions, and compare numerical and exact values at a fixed physical time. To measure spatial order, halve $h$ while making temporal error negligible; the <error norm> should fall by a factor of four. To measure temporal order, halve $k$ on a sufficiently fine fixed spatial grid and again expect a factor of four. Report observed order $p=\log_2(E_h/E_{h/2})$ and use several refinements to identify the asymptotic regime.

Check more than one smooth mode, boundary implementation and expected conservation or decay properties. A <manufactured solution> with an added known source can test nontrivial boundaries and forcing. Repeat in maximum and root-mean-square <norms>, stop before finite-precision or solver tolerances dominate, and separate stability failures from errors in convergence order. A single small error or a single visually smooth solution does not establish accuracy.