We present numerical methods for analysis of the stability of two-dimensional steady viscoelastic flows to small-amplitude, two-dimensional and three-dimensional disturbances based on finite element calculations of the steady base flow and the perturbation. Direct time integration of the linearized equations of motion and iterative calculation of the most dangerous components of the eigenspectrum are tested. Finite element discretizations based on the DEVSS-G finite element discretization with Newton's method used to compute steady-state solutions. Two different rime integration schemes are tested for computing the time evolution of general, random disturbances: a theta-method operator-splitting scheme and a fourth-order Runge-Kutta method. For both time integrators, time stepping is decoupled into a solution of a modified Stokes problem and an evaluation of the time-dependent constitutive equation. The overall efficiency of both methods is extremely high, as is the potential for implementation on parallel computers. An algorithm also is presented for calculating eigenvalues with the largest real part that combines time integration of the Linearized equations with a Krylov subspace method to accelerate the calculation of the eigenvalues. Although this method does not dramatically reduce the computational cost over use of time integration alone, it does provide a more complete analysis of the eigenspectrum. For both direct time integration and the hybrid time integration/Krylov calculation, the stability results for cylindrical Couette flow show quantitative agreement with the eigenvalues calculated by using other methods of analysis [M. Avgousti, A.N. Beris, J. Non-Newtonian Fluid Mech. 50 (1993) 225-251]. Contrary to the results in our previous paper [R. Sureshkumar, M.D. Smith, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 82 (1999) 57-104], we find that the flow of an Oldroyd-B fluid through a closely-spaced cylinder array is stable to two-dimensional perturbations. However, allowing the perturbations to be three-dimensional and considering an isolated cylinder does not alter the conclusions of our earlier study [R. Sureshkumar, M.D. Smith, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 82 (1999) 57-104] of the two-dimensional stability of widely-spaced arrays of cylinders; the flow around an isolated cylinder is computed to be stable for all values of the Weissenberg number obtainable with these calculations, We less than or equal to 0.75.