Discretization of the Stokes equations produces a symmetric indefinite system of linear equations. For stable discretizations, a variety of numerical methods have been proposed that have rates of convergence independent of the mesh size used in the discretization. In this paper, we compare the performance of four such methods: variants of the Uzawa, preconditioned conjugate gradient, preconditioned conjugate residual, and multigrid methods, for solving several two-dimensional model problems. The results indicate that where it is applicable, multigrid with smoothing based on incomplete factorizaton is more efficient than the other methods, but typically by no more than a factor of two. The conjugate residual method has the advantages of being both independent of iteration parameters and widely applicable. (Also cross-referenced as UMIACS-TR-94-76)