In this paper, numerical modelling of transient viscoelastic flows in the context of a finite volume formulation is discussed with focus on numerical accuracy, efficiency and stability of the iterative solution algorithms used in solving the coupled system of the discretized governing equations. Generalized three-level schemes are used for temporal discretization, which give at least second-order temporal discretization accuracy for a convection-diffusion equation. However, we will show that numerical stability and efficiency of calculations for the coupled system of the governing equations can be poor if an improper iterative solution algorithm is used. It is found that the explicit solution algorithms in which the parent equations are discretized explicitly in time is not suitable for transient viscoelastic calculations due to the possible spatial oscillations whenever the order of spatial discretization for hyperbolic constitutive equations is higher than the first-order. On the other hand, with the implicit solution algorithms for truly transient flow calculations, only when the indicative errors in the iterative solution process are at a compatible (the same or smaller) order with that embodied in the discretization process, can we take the temporal discretization order to be the accuracy of the final solutions for the coupled system. In this sense, the SIMPLEST and PISO algorithms are two of the most efficient and stable algorithms available. Our analysis has been confirmed by numerical calculations carried out for the start-up and decay of Poiseuille flows between two parallel plates for Newtonian fluid and for viscoelastic fluids with a solvent viscosity described by the Oldroyd-B model. We demonstate that analytical solutions can be reproduced smoothly without limitation on Weissenberg number. Finally, we analyze and numerically demonstrate that with the decoupled method, for accurate simulations of the transient viscoelastic flow with inertia, the governing equations have to be decoupled and solved strictly according to the type of each equation. The stabilizing measures developed for steady state modelling, such as 'both sides diffusion' (BSD), elastic viscous stress split (EVSS), and adaptive viscous split stress (AVSS) may lead to change of type in the equation from a purely hyperbolic one to a parabolic one for a purely elastic (Maxwell) fluid, or lead to overdiffusion in transient velocity field calculations for the viscoelastic fluid with a solvent viscosity. Therefore, these kinds of measures should not be taken, and the 'viscous formulation' is more proper. (C) 2004 Elsevier B.V. All rights reserved.