A new method for stabilizing viscoelastic flows is proposed suitable for high-order discretizations. It employs a mode-dependent diffusion operator that guarantees monotonicity while maintaining the formal accuracy of the discretization. Other features of the method are: a high-order time-splitting scheme, modal spectral element expansions on a single grid, and the use of a finitely extensible non-linear elastic-Peterlin (FENE-P) model. The converuence of the method is established through analytic examples and benchmark problems in two and three dimensions, and unsteady flow past a three-dimensional (313) ellipsoid is studied at high Reynolds number. (C) 2003 Elsevier B.V. All rights reserved.