  Accretion disks are three-dimensional, turbulent, often self-gravitating,
magnetohydrodynamic flows, which can be modeled in detail with numerical
simulations. In this paper, we present a new algorithm that is based on a
spectral decomposition method to simulate such flows. Because of the high order
of the method, we can solve the induction equation in terms of the magnetic
potential and, therefore, ensure trivially that the magnetic fields in the
numerical solution are divergence free. The spectral method also suffers
minimally from numerical dissipation and allows for an easy implementation of
models for sub-grid physics. Both properties make our method ideal for studying
MHD turbulent flows such as those found in accretion disks around compact
objects. We verify our algorithm with a series of standard tests and use it to
show the development of MHD turbulnce in a simulation of an accretion disk.
Finally, we study the evolution and saturation of the power spectrum of MHD
turbulence driven by the magnetorotational instability.
