A numerical solution of a partial differential equation (PDE) with expected large solution gradients demands a high level of resolution in approximating temporal and spatial derivatives. Low-order methods, first- and second-order ones (reviewed in Appendix A ), provide limited resolution for a given grid spacing since their errors scale as \(\mathcal {O}(\tau )\) or \(\mathcal {O}(\tau ^2)\) for a temporal scheme and \(\mathcal {O}(h)\) or \(\mathcal {O}(h^2)\) for a spatial one. In other words, the rate of convergence of the first- and second-order methods is low; hence, high-fidelity computations require an exceedingly high level of grid refinements. The high number of grid points in space and time can lead to prohibitively expensive computations.