An efficient boundary element solution of the motion of inelastic non-Newtonian fluids at low Reynolds number is presented in this paper. For the numerical solution all the domain integrals of the boundary element formulation have been transformed into equivalent boundary integrals by means of the dual reciprocity method (DRM). To achieve an accurate approximation of the non-linear and non-Newtonian terms two major improvements have been made to the DRM, namely the use of augmented thin plate splines as interpolation functions, and the partition of the entire domain into smaller subregions or domain decomposition. In each subregion or domain element the DRM was applied together with some additional equations that ensure continuity on the interfaces between adjacent subdomains. After applying the boundary conditions the final systems of equations will be sparse and the approximation of the nonlinear terms will be more localised than in the traditional DRM. This new method known as multidomain dual reciprocity (MD-DRM) has been used to solve several non-Newtonian problems including the pressure driven flow of a power law fluid, the Couette flow and two simulations of industrial polymer mixers.