A finite volume method for solving Navier-Stokes equations in cases of laminar incompressible flows was developed. The method is based on a new-proposed pressure-velocity coupling algorithm. The fundamental difference of the new algorithm, compared to the most frequently used SIMPLE algorithm, is that in the new algorithm, the continuity equation is satisfied in each iteration. Within the new algorithm, the final solution is approached by correcting the cell face mass flux until the pressure gradient field becomes irrotational. The new algorithm (implemented in a in-house computer code) and the SIMPLE algorithm (implemented in the commercial software package, FLUENT) are compared in terms of convergence rate in few typical cases of laminar flow ranging from low to high Reynolds number, and in cases of natural convection covering the range from low to high Rayleigh number. In accordance with the obtained results, it is concluded that in the cases of inviscid and viscous flow: (1) The convergence rate of the new algorithm is significantly higher, and the user need not to prescribe pressure and velocity under-relaxation factors on which this rate would depend, as is the case with the SIMPLE algorithm; (2) The convergence rate of the new algorithm is practically independent on mesh density (mesh size), in contrast to the convergence rate of the SIMPLE algorithm, which decreases with increasing mesh density. Rate of convergence of the new algorithm slightly decreases with increasing complexity of flow pattern (the appearance of vortices and zones of flow separation) and with applying a higher order differencing scheme (due to deferred correction approach). Here, when applying the higher order scheme, the convergence rate increases with increasing mesh density; (3) The new algorithm shows better efficiency than SIMPLE, regardless to possibly suboptimal coding of the new algorithm when compared to the SIMPLE algorithm, professionally coded within the FLUENT software package. In cases of natural convection the number of iterations required to achieve the prescribed solution accuracy is greater with respect to the cases of pure inviscid and laminar flows, due to partially explicit treatment of the buoyancy force. Here, in order to increase the convergence rate of the new algorithm, the equations for the mass flux and temperature corrections should be treated in a coupled manner. The number of iterations required in new algorithm increases with the Rayleigh number, but here also, throughout the whole range of Rayleigh number, the new algorithm shows better performance than the SIMPLE algorithm. The new algorithm results in a linear system with less unknowns, but with a denser matrix than in the case of SIMPLE algorithm, thus it is more memory demanding. The drawback of the new algorithm is that this matrix does not possess the properties of the well known M-matrix, as in the case of the SIMPLE algorithm, which requires additional research in finding of an efficient iterative solver for such type of matrix.