The numerical simulations of non-Newtonian flows become difficult as the Weissenberg number increases. The main objective of this study is to generate a robust and efficient matrix-free Inexact Newton-Krylov solver (IN-GMRES) which can deal with the high nonlinearities arising from the increased Weissenberg numbers. In order to achieve this aim, non-Newtonian flows of rotating surface driven cylindrical enclosure problem is numerically investigated by using three differential viscoelastic constitutive relations namely Upper Convected Maxwell (UCM), Oldroyd B and Giesekus. The results obtained by IN-GMRES solver are validated and compared with the POLYFLOW simulations. Additionally, the selection of the constitutive relation, effects of the Weissenberg number and effects of the Reynolds number are studied. The simulations indicate that the generated algorithm is capable of solving higher Weissenberg number problems (up to the Elasticity number limit of 130) when compared to the previous studies. Furthermore, it is shown that with the increasing Weissenberg number, the reversed flow can be observed in the flow domain and in some cases, depending on the Reynolds number, re-formation of the Newtonian like flow is possible at high Weissenberg numbers.