The analysis of steady incompressible flow is based on the equation of motion and the continuity constraint together with appropriate boundary conditions. Besides the inertia term, also the extra stress term becomes non-linear in the case of generalized Newtonian fluids. By application of mixed finite-element methods the conservation equations are reduced to a system of non-linear algebraic equations. At each step of a Newton-Raphson scheme a linear subproblem with a non-symmetric and indefinite coefficient matrix has to be solved. A multigrid method for these systems is presented, where the smoothing iteration is a modification of Uzawa's algorithm based on an augmented Lagrangian functional. The method is applied to the numerical simulation of vortex breakdown in a cylindrical vessel, where the computations are performed for real shear-thinning fluids. Dimensionless groups, a Reynolds number and a non-Newtonian parameter are introduced, and the dependence of the convergence rate on these parameters is discussed.