To improve image quality in susceptibility-weighted MR imaging, it is important to correct for the effects of field inhomogeneity. In particular, susceptibility differences between air and tissue induce magnetic field nonuniformity; often those susceptibility effects have nonzero through-plane gradients that lead to spin dephasing across the slice within each voxel. If uncorrected, these through-plane gradients cause signal loss in the reconstructed images. Several methods exist for reconstructing MR images with compensation for field inhomogeneity, but most of these methods, even the model-based iterative ones, treat the inhomogeneity within each voxel as being a constant, thus ignoring the through-plane gradient effects. This paper describes a model-based iterative method for reconstructing MR images with compensation for field inhomogeneity that accounts for the slice profile and the throughplane gradients of the field inhomogeneity (assumed to be determined by a pre-scan). In particular, this paper describes an accelerated algorithm for implementing the forward model and its adjoint as needed in a conjugate gradient algorithm for iterative MR image reconstruction.