  Time dependent signals in experimental techniques such as Nuclear Magnetic
Resonance (NMR) and Muon Spin Relaxation (muSR) are often the result of an
ensemble average over many microscopical dynamical processes. While there are a
number of functions used to fit these signals, they are often valid only in
specific regimes, and almost never properly describe the "spectral diffusion"
regime, in which the dynamics happen on time scales comparable to the
characteristic frequencies of the system. Full treatment of these problems
would require one to carry out a path integral over all possible realizations
of the dynamics of the time dependent Hamiltonian.
  In this paper we present a numerical approach that can potentially be used to
solve such time evolution problems, and we benchmark it against a Monte Carlo
simulations of the same problems. The approach can be used for any sort of
dynamics, but is especially powerful for any dynamics that can be approximated
as Markov processes, in which the dynamics at each step only depend on the
previous state of the system. The approach is used to average both classical
and quantum observables; in the latter case, a formalism making use of
Liouvillians and density matrices is used.
