Summary: We study the application of Newton-Krylov methods for the steady-state solution of the tokamak edge plasma fluid equations. This highly nonlinear system of two-dimensional convection-diffusion-reaction partial differential equations describes the boundary layer of a tokamak fusion reactor. These equations are characterized by multiple time and spatial scales. We use Newton’s method to linearize the nonlinear system of equations resulting from the implicit, finite volume discretization of the governing partial differential equations. The resulting linear systems are neither symmetric nor positive definite, and are poorly conditioned. A variety of preconditioned Krylov iterative techniques are employed to solve these linear systems, and we investigate both standard and matrix-free implementations. A number of pseudo-transient continuation methods are investigated to increase the radius of convergence. While this system of equations describes a specific application, the general algorithm should benefit other applications requiring the solution of general reacting flow type equations.