  In this paper we revisit an approach pioneered by Auchmuty to approximate
solutions of the Laplace- Robin boundary value problem. We demonstrate the
efficacy of this approach on a large class of non-tensorial domains, in
contrast with other spectral approaches for such problems.
  We establish a spectral approximation theorem showing an exponential fast
numerical evaluation with regards to the number of Steklov eigenfunctions used,
for smooth domains and smooth boundary data. A polynomial fast numerical
evaluation is observed for either non-smooth domains or non-smooth boundary
data. We additionally prove a new result on the regularity of the Steklov
eigenfunctions, depending on the regularity of the domain boundary.
  We describe three numerical methods to compute Steklov eigenfunctions.
