A front-tracking method is developed for direct numerical simulations of viscoelastic two-phase systems in which one or both phases could be viscoelastic. One set of governing equations is written for the whole computational domain and different phases are treated as a single fluid with variable material and rheological properties. The interface is tracked explicitly using a Lagrangian grid while the flow equations are solved on a fixed Eulerian grid. The surface tension is computed at the interface using the Lagrangian grid and included into the momentum equations as a body force. The Oldroyd-B, FENE-CR and FENE-MCR models are employed to model the viscoelasticity. The viscoelastic model equations are solved fully coupled with the flow equations within the front-tracking framework. A fifth-order WENO scheme is used to approximate the convective terms in the viscoelastic model equations and second-order central differences are used for all other spatial derivatives. A log-conformation method-is employed to alleviate the high Weissenberg number problem (HWNP) and found to be stable and very robust for a wide range of Weissenberg numbers. The method has been first validated for various benchmark single-phase and two-phase viscoelastic flow problems. Then it has been applied to study motion and deformation of viscoelastic two-phase systems in a pressure-driven flow through a capillary tube with a sudden contraction and expansion. The method has been demonstrated to be grid convergent with second-order spatial accuracy for all the cases considered in this paper. (C) 2015 Elsevier B.V. All rights reserved.