In this thesis, a new finite element method, 4-field Galerkin/Least-Squares method, is presented to solve viscoelastic flow problems. The 4-field GLS naturally includes the SUPG and PSPG terms to stabilize the oscillations caused by advection-dominated terms. In addition, it introduces a new variable L = ∇v, so that the second order derivative of v is avoided, and the basis functions can be chosen as piecewise linear functions. This feature substantially enlarges the space of the basis and weighting functions. The Galerkin terms in this formulation guarantee that the traction term n·T appears naturally by integration by part, which serves as an important boundary condition for free surface flow. Moreover, the 4-field GLS successfully circumvents the LBB condition on velocity and conformation fields. The 4-field GLS is tested with a carefully defined set of benchmark problems for both Newtonian and non-Newtonian fluid. It is found to be robust, accurate and efficient.