This paper presents a boundary-element-only method (BEM) for the calculation of generalised Newtonian fluid (GNF) flows. The volume integral arising from non-linear effects is approximated via a particular solution technique. Multilayer perceptron networks (MLPN) and radial basis function networks (RBFN) are used for global approximation of field variables and hence volume discretisation is not required. The iterative numerical formulation is achieved by viewing the material as being composed of a Newtonian base (artificially assigned with a constant, but maybe different from subdomain to subdomain, viscosity) and the remaining component, which is accordingly defined from the original constitutive equation. This decoupling of the non-linear effects allows a Picard-type iterative procedure to be employed by treating the non-linear term as a known forcing function. However, convergence is sensitive to the estimate of this forcing function and an adaptive subregioning of the domain is adopted to control the accuracy of the estimate of this non-linear term. The criterion for subregioning is that the velocity gradient should not vary significantly in each subdomain. This strategy enables convergence of the present method (BEM-NN) at power-law index as low as 0.2 for the difficult power law fluid. The use of MLPNs (instead of single layer perceptrons) and RBFNs is another contributing factor to the improved convergence performance. The overall scheme is very suitable for coarse-grain parallelisation as each subdomain can be independently analysed within an iteration. Furthermore, within each subdomain process, there are other parallelisable computations. The present method is verified with circular Couette and planar Poiseuille flows of the power-law, Carreau–Yasuda and Cross fluids.