To increase the accuracy of computations by the method of collocations and least squares (CLS) a generalization of this method is proposed for the case of a non-uniform logically rectangular grid. The main work formulas of the CLS method on non-uniform grid, including the formulas implementing the prolongation operator on a non-uniform grid at the use of a multigrid complex are obtained with the aid of the computer algebra system (CAS) Mathematica. The proposed method has been applied for the numerical solution of two-dimensional stationary Navier–Stokes equations governing the laminar flows of viscous incompressible fluids. On a smooth test solution, the application of a non-uniform grid has enabled a 47-fold reduction of the solution error in comparison with the uniform grid case. At the solution of the problem involving singularities – the lid-driven cavity flow – the error of the solution obtained by the CLS method was reduced by the factors from 2.65 to 3.05 depending on the Reynolds number value.