A smoothed particle hydrodynamics approach is utilized to model a non-Newtonian fluid with a spatially varying viscosity. In the limit of constant viscosity, this approach recovers an earlier model for Newtonian fluids of Espaol and Revenga (Phys Rev E 67:026705, 2003). Results are compared with numerical solutions of the general Navier-Strokes equation using the "regularized" Bingham model of Papanastasiou (J Rheol 31:385-404, 1987) that has a shear-rate-dependent viscosity. As an application of this model, the effect of having a non-Newtonian fluid matrix, with a shear-rate-dependent viscosity in a moderately dense suspension, is examined. Simulation results are then compared with experiments on mono-size silica spheres in a shear-thinning fluid and for sand in a calcium carbonate paste. Excellent agreement is found between simulation and experiment. These results indicate that measurements of the shear viscosity of simple shear-rate-dependent non-Newtonian fluids may be used in simulation to predict the viscosity of concentrated suspensions having the same matrix fluid.