We consider the bidimensional Stokes problem for incompressible fluids in stream function-vorticity formulation. For this problem, the classical finite elements method of degree one converges only in 𝒪(h) for the quadratic norm of the vorticity, if the domain is convex and the solution regular. We propose to use harmonic functions obtained by a simple layer potential to approach vorticity along the boundary. Numerical results are very satisfying and we prove that this new numerical scheme leads to an error of order 𝒪(h) for the natural norm of the vorticity and under more regularity assumptions from 𝒪(h 3/2 ) to 𝒪(h 2 ) for the quadratic norm of the vorticity.