The well known Reid-Harris expansions, applied to the stream function formulation, and the Projection-Diffusion Chebyshev Stokes solver, in primitive variables, are used to compute the fundamental Stokes eigenmodes of each of the symmetry families characterizing the Stokes solutions in the square. The numerical accuracy of both methods, applied with several discretizations, are compared, for both the eigenvalues and the main features of the corresponding eigenmodes. The Chebyshev approach is by far the most efficient, even though the associated solver does not provide a divergence free velocity but asymptotically.