  Single-stage or single-step high-order temporal discretizations of partial
differential equations (PDEs) have shown great promise in delivering high-order
accuracy in time with efficient use of computational resources. There has been
much success in developing such methods for finite volume method (FVM)
discretizations of PDEs. The Picard Integral formulation (PIF) has recently
made such single-stage temporal methods accessible for finite difference method
(FDM) discretizations. PIF methods rely on the so-called Lax-Wendroff
procedures to tightly couple spatial and temporal derivatives through the
governing PDE system to construct high-order Taylor series expansions in time.
Going to higher than third order in time requires the calculation of
Jacobian-like derivative tensor-vector contractions of an increasingly larger
degree, greatly adding to the complexity of such schemes. To that end, we
present in this paper a method for calculating these tensor contractions
through a recursive application of a discrete Jacobian operator that readily
and efficiently computes the needed contractions entirely agnostic of the
system of partial differential equations (PDEs) being solved.
