This paper addresses the analysis of spectrum and pseudospectrum of the linearized Navier-Stokes operator from the numerical point of view. The pseudospectrum plays a crucial role in linear hydrodynamic stability theory and is closely related to the non-normality of the underlying differential operator and the matrices resulting from its discretization. This concept offers an explanation for experimentally observed instability in situations when eigenvalue-based linear stability analysis would predict stability. Hence the reliable numerical computation of the pseudospectrum is of practical importance particularly in situations when the stationary "base flow" is not analytically but only computationally given. The proposed algorithm is based on a finite element discretization of the continuous eigenvalue problem and uses an Arnoldi-type method involving a multigrid component. Its performance is investigated theoretically as well as practically at several two-dimensional test examples such as the linearized Burgers equations and various problems governed by the Navier-Stokes equations for incompressible flow.