A method for computing large numbers of eigenvalues of self-adjoint elliptic operators [J. Comput. Phys. 184, 321 (2003)] is tested in numerical studies of the spectral properties of quantum billiards. To this extent, we study a time-reversal invariant quantum billiard of threefold symmetry, that undergoes a transformation in its symmetry properties from C(3v) to C3 . Thereby a transition from Gaussian orthogonal to Gaussian unitary ensemble statistics is observed, verifying earlier experimental indications and theoretical predictions. At the same time our numerical ansatz is shown to be applicable to arbitrary billiard shapes.