Abstract. This paper develops a nested iteration algorithm to solve time-dependent nonlinear systems of partial differential equations. For each time step, Newton's method is used to form approximate solutions from a sequence of nested spaces, where the resolution of the approximations increases as the algorithm progresses. Nested iteration results in most of the iterations being performed on coarser grids, where minimal work is needed to reduce error to the level of discretization error. The approximate solution on a given coarse grid is interpolated to a refined grid and used as an initial guess for the problem posed there. The approximation is then already close enough to the solution on the current grid that a minimal amount of work is needed to solve the refined problem, due to the rapid convergence of Newton's method near a solution. The paper develops an algorithm that attempts to optimize accuracy-percomputational-cost on each grid, so that essentially no unnecessary work is done on any grid. The nested iteration algorithm is then applied to a reduced 2D model of the incompressible, resistive magnetohydrodynamic (MHD) equations. Using this algorithm on the MHD equations in the context of a first-order system least-squares (FOSLS) finite element discretization and algebraic multigrid (AMG) to solve the linearized systems, instabilities in a model tokamak fusion reactor are simulated. Numerical results show that this highly complex nonlinear problem is solved in an equivalent of 30-80 fine-grid relaxation sweeps per time step.Key words. magnetohydrodynamics, algebraic multigrid, FOSLS, nested iteration AMS subject classifications. 65F10, 65M55, 76W05Acknowledgments. This work was sponsored by the Department of Energy under grant numbers DE-FG02-03ER25574 and DE-FC02-06ER25784, Lawrence Livermore National Laboratory under contract numbers B568677, and the National Science Foundation under grant numbers DMS-0621199, DMS-0749317, and DMS-0811275.1. Introduction. This paper applies a nested iteration algorithm to a set of incompressible resistive Magnetohydrodynamic (MHD) equations. A first-order systems least-squares (FOSLS) [11,12] finite element discretization is used. The main focus of this paper is to show that nested iteration is a crucial component to solving complicated systems of nonlinear partial differential equations. As is shown for the MHD equations, nested iteration allows for a complicated system of nonlinear equations to be solved to within discretization accuracy in only a handful of work units (WUs). A WU is defined as the amount of computation required to perform one fine-grid relaxation sweep. The main idea is that most of the work is done on coarse grids where computation is much cheaper. The approximation to the solution on these coarse grids is interpolated up to successively finer grids, where linearizations and algebraic solves are applied. The process is then continued to finer grids until a desired error tolerance or resolution is reached. The result is that coarse grids give fairly...