  As a follow up to \cite{Causley2013}, we provide a detailed description of
the numerical implementation of an O(N), A-stable, second order accurate
solution of the wave equation, constructed from semi-discrete boundary value
problems. We improve on the previous algorithm by replacing the Lax-type
correction used in \cite{Causley2013}, which was necessary for convergence when
$\Delta t < \Delta x/c$, with a more accurate spatial quadrature, which we
prove is convergent.
  We also demonstrate that the resulting solver remains fast even in the case
of unstructured meshes, can incorporate domain decomposition, and allows for
the implementation of Dirichlet, Neumann, periodic and outflow boundary
conditions.
  Building upon results for the 1d formulation, we utilize alternate direction
implicit (ADI) splitting to achieve a fast O(N) solver in higher spatial
dimensions. Our solver is built upon line objects and, combined with the
flexibility of the integral solver, allows us to solve problems on arbitrary
spatial domains, by embedding the boundary in a regular Cartesian mesh. Our
solver is designed to couple with particle codes, where scale separation is an
issue. We therefore demonstrate the ability of our solver to take time steps
well beyond that of the Courant-Friedrichs-Lewy (CFL) stability limit of
explicit codes.
