Existing methods for surface quadrangulation cannot ensure accurate alignment with feature or boundary curves and tight control of local element size, which are important requirements in many numerical applications (e.g., FEA). Some methods rely on a prescribed direction field to guide quadrangulation for feature alignment, but such a direction field may conflict with a desired density field, thus making it difficult to control the element size. We propose a new spectral method that achieves both accurate feature curve alignment and tight control of local element size according to aThe work of W.given density field. Specifically, the following three technical contributions are made. First, to make the quadrangulation align accurately with feature curves or surface boundary curves, we introduce novel boundary conditions for wave-like functions that satisfy the Helmholtz equation approximately in the least squares sense. Such functions, called quasi-eigenfunctions, are computed efficiently as the solutions to a variational problem. Second, the mesh element size is effectively controlled by locally modulating the Laplace operator in the Helmholtz equation according to a given density field. Third, to improve robustness, we propose a novel scheme to minimize the vibration difference of the quasi-eigenfunction in two orthogonal directions. It is demonstrated by extensive experiments that our method outperforms previous methods in generating feature-aligned quadrilateral meshes with tight control of local elememt size. We further present some preliminary results to show that our method can be extended to generating hex-dominant volume meshes.