A new mixed finite element method is presented that has improved numerical stability and has similar numerical accuracy compared to the EVSS/FEM developed by Rajagopalan et al. (J. Non-Newtonian Fluid Mech. 36 (1990) 159-192). The new method, denoted EVSS-G/ FEM, is based on separate bilinear interpolation of all relevant components of the velocity gradient tenser. Lagrangian bilinear approximations of this variable are formed by least-squares interpolation of the derivatives of a Lagrangian biquadratic approximation to the velocity field that is generated by solution of the momentum and continuity equations using a standard mixed formulation for velocity and pressure. Elastic components of the stress are computed by solving the constitutive equation using either the streamline upwinding (SU) or streamline upwind Petrov-Galerkin (SUPG) methods with bilinear interpolation. The motivation for the interpolation of the velocity gradients and for the stress interpolation is to force the compatibility between these variables that is needed to satisfy differential viscoelastic constitutive equations in the limit where the velocity field vanishes, as it does near solid boundaries. The enhanced numerical stability of the EVSS-G/FEM over the EVSS/FEM formulation is demonstrated for the calculation of the linear stability of planar Couette flow of a UCM fluid, where closed form results exist for the linear stability of the steady-state flow; hence, any instability detected in the numerical solution is a result of either the finite element approximation or the time integration method. The second possibility is eliminated by using a fully implicit time integration method. Calculations with the EVSS-G/FEM are stable for De > 100, whereas calculations with the EVSS/FEM become numerically unstable for De > 5. The EVSS-G formulation is combined with finite element discretizations that incorporate local adaptive mesh refinement for increasing the accuracy of the calculations in regions where the stress and velocity gradients change rapidly. The accuracy of the EVSS-G/FEM is demonstrated by steady-state calculations for the flow between eccentric rotating cylinders, flow through a wavy-walled tube, and for flow through a square array of cylinders; each has been used before as a test problem for the accuracy of viscoelastic flow calculations. As in previous calculations, discretization of the constitutive equation by SUPG gives superior accuracy at high values of De compared to the application of SU.