The solution of the equations governing the steady incompressible slow viscous fluid flow is analysed using a novel technique based on a Laplacian decomposition instead of the more traditional approaches based on the biharmonic streamfunction formulation or the velocity-pressure formulation. This results in the need to solve the Laplace equations for the pressure and other auxiliary harmonic functions which arise from the ideas of Almansi's decomposition. These equations, which become coupled through the boundary conditions, are numerically solved using the boundary element method (BEM). Results both on the boundary and inside the solution domain are presented and discussed for a simple benchmark test example and a few applications in smooth and non-smooth geometries in order to illustrate that the Laplacian decomposition in combination with BEM provides an efficient technique, in terms of accuracy and convergence, to investigate numerically a Stokes flow.