In this thesis, we develop an adaptive mesh refinement (AMR) code including magnetic fields, and use it to perform high resolution simulations of magnetized molecular clouds. The purpose of these simulations is to study present day star formation in the presence of turbulence and magnetic fields. We first present MHDEnzo, the extension of the cosmology and astrophysics code Enzo to include the effects magnetic fields. We use a higher order Godunov Riemann solver for the computation of interface fluxes; constrained transport to compute the electric field from those interface fluxes, which advances the induction equation in a divergence free manner; divergence free reconstruction technique to interpolate the magnetic fields to fine grids; operator splitting to include gravity and cosmological expansion. We present a series of test problems to demonstrate the quality of solution achieved. Additionally, we present several other solvers that were developed along the way. Finally we present the results from several AMR simulations that study isothermal turbulence in the presence of magnetic fields and self gravity. Ten simulations with initial Mach number 8.9 were studied varying several parameters; virial parameter $\alpha$ from 0.52 to 3.1; whether they were continuously stirred or allowed to decay ; and the number of refinement levels (4 or 6). Measurements of the density probability density function (PDF) were made, showing both the expected log normal distribution and an additional power law. Measurements of the line of sight magnetic field vs. column density are done, giving excellent agreement with recent observations. The line width vs. size relationship is measured and compared with good agreement to observations, reproducing both turbulent and collapse signatures. The core mass distribution is measured and agrees well with observations of Serpens and Perseus core samples, but the power-law distribution in Ophiuchus is not reproduced by our simulations. Finally we attempt to make contact with recent theoretical predictions of the star formation rate. Our measured rate is significantly higher than predicted, indicating that our root grid resolution is likely too low. Nonetheless, the simulations presented here are the first of their kind, and the general agreement with observations indicates the promise of our approach. We conclude by outlining future work which will explore numerical systematics more fully, and make detailed contact with observations.