This paper reports the developments made to improve the numerical stability of the open-source finite-volume computational library OpenFOAM (R) developed for the numerical computation of viscoelastic fluid flows described by differential constitutive models. The improvements are based on the modification of the both-sides diffusion technique, named improved both-sides diffusion (iBSD), which promotes the coupling between velocity and stress fields. Calculations for two benchmark 2D case studies of an upper-convected Maxwell (UCM) fluid are presented and compared with literature results, namely the 4:1 planar contraction flow and the flow around a confined cylinder. The results obtained for the first case are computed in five meshes with different refinement levels and are compared with literature results. In this case study it was possible to achieve steady-state converged solutions in the range of Deborah numbers tested, De = {0, 1, 2, 3, 4, 5}, for all meshes. The corner vortex size predictions agree well with the literature and a relative error below 0.6% is obtained for De <= 5. In the flow around a confined cylinder, steady-state converged solutions were obtained in the range of Deborah numbers tested, De = {0, 0.3, 0.6, 0.8), in four consecutively refined meshes. The predictions of the drag coefficient on the cylinder are similar to reference data with a relative error below 0.08%. For both test cases the developed numerical method was shown to have a convergence order between 1 and 2, in general very close to the latter. Moreover, the results presented for both case studies clearly extend the previous ones available in the literature in terms of accuracy. This was a direct consequence of the capability of performing the calculation with more refined meshes, than the ones employed before.