This paper is concerned with the accurate numerical approximation of the spectral properties of the biharmonic operator on various domains in two dimensions. A number of analytic results concerning the eigenfunctions of this operator are summarized and their implications for numerical approximation are discussed. In particular, the asymptotic behaviour of the first eigenfunction is studied since it is known that this has an unbounded number of oscillations when approaching certain types of corners on domain boundaries. Recent computational results of Bjørstad & Tjøstheim, using a highly accurate spectral Legendre–Galerkin method, have demonstrated that a number of these sign changes may be accurately computed on a square domain provided sufficient care is taken with the numerical method. We demonstrate that similar accuracy is also achieved using an unstructured finite–element solver which may be applied to problems on domains with arbitrary geometries. A number of results obtained from this mixed finite–element approach are then presented for a variety of domains. These include a family of circular sector regions, for which the oscillatory behaviour is studied as a function of the internal angle, and another family of (symmetric and non–convex) domains, for which the parity of the least eigenfunction is investigated. The paper not only verifies existing asymptotic theory, but also allows us to make a new conjecture concerning the eigenfunctions of the biharmonic operator.