In the framework of finite element analysis we propose fast and robust time integration scheme for viscoelastic fluid (the Oldroyd-B and Leonov models) flow by way of efficient decoupling of equations. Developed algorithms of the $1^{st}$ and $2^{nd}$ order are shown to disclose convergence characteristics equivalent to conventional methods of corresponding order when applied to 1D poiseuille and 2D creeping contraction flow problems. In comparison with fully coupled implicit technique, they notably enhance the computation speed. For the time dependent flow modeling with pressure difference imposed slightly below the steady limit, current as well as conventional approximation scheme has demonstrated fluctuating solution without approaching the steady state. From the result, we may conclude that the existence of upper limit for convergent steady solution implies flow transition to highly elastic time-fluctuating field without steady asymptotic. It is presumably associated with some real unstable elastic flow like re-entrant vortex oscillation and extrudate distortion outside the channel outlet.