A robust and efficient numerical scheme has been developed for the solution of the finite-differenced pressure linked fluid flow equations. The algorithm solves the set of nonlinear simultaneous equations by a combination of Newton's method and efficient sparse matrix techniques. In tests on typical recirculating flows the method is rapidly convergent. The method does not require any under-relaxation or other convergence-enhancing techniques employed in iterative schemes. It is currently described for two-dimensional steady state flows but is extendible to three dimensions and mildly time-varying flows. The method is robust to changes in Reynolds number, grid aspect ratio, and mesh size. This paper reports the algorithm and the results of calculations performed.