An incremental difference formulation (IDF) has been developed to further improve the efficiency and stability of the operator splitting algorithm for viscoelastic flows, taking advantage of the pseudo linearity of the governing equations. Using this new IDF/OS scheme and up to 200000 degrees of freedom, high resolution mesh-convergent solutions for the falling sphere problem are obtained up to Wi = 2.8, with a minimum value of K = 4.02 predicted at Wi = 2.2. Beyond Wi = 2.8, the solutions become significantly mesh-dependent, although steady convergence with individual meshes are still possible. Well structured quadrilateral elements lead to better accuracy than the entirely unstructured triangle meshes previously used. There is some local mesh-dependence of T-zz on the centreline starting at a finite distance downstream from the rear stagnation point, even for globally mesh-convergent solutions. This local mesh-dependence is most likely due to the presence of a thin boundary layer forming along the centreline in the stress wake, where T-zz exhibits a sharp radial gradient and its local maximum is not on the centreline but slightly off the centreline inside the boundary layer. Further mesh refinement in the radial direction, perhaps reducing the size by an order of magnitude, is needed to study the nature of this boundary layer.