  In several recent works \cite{Causley2013a}, \cite{Causley2013}, we developed
a new second order, A-stable approach to wave propagation problems based on the
method of lines transpose (MOL$^T$) formulation combined with alternating
direction implicit (ADI) schemes. In this work, we present several important
modifications to our work, and thus obtain a family of wave solvers which are
unconditionally stable, accurate of order 2P, and require $O(P^d N)$ operations
per time step, where $N$ is the number of spatial points, and $d$ the number of
spatial dimensions. We obtain these schemes by including higher derivatives of
the solution, rather than increasing the number of time levels. The novel
aspect of our approach is that the higher derivatives are constructed using
successive applications of the convolution operator.
  We develop these schemes in one spatial dimension, and then extend the
results to higher dimensions, by reformulating the ADI scheme to include
recursive convolution. Thus, we retain a fast, unconditionally stable scheme,
which does not suffer from the large dispersion errors characteristic to the
ADI method. We demonstrate the utility of the method by applying it to a host
of wave propagation problems. This method holds great promise for developing
higher order, parallelizable algorithms for solving hyperbolic PDEs, and can
also be extended to parabolic PDEs.
