Taylor bubbles are a feature of the slug flow regime in gas–liquid flows in vertical pipes. Their dynamics exhibits a number of transitions such as symmetry breaking in the bubble shape and wake when rising in downward flowing and stagnant liquids, respectively, as well as breakup in sufficiently turbulent environments. Motivated by the need to examine the stability of a Taylor bubble in liquids, a systematic numerical study of a steadily moving Taylor bubble in stagnant and flowing liquids is carried out, characterised by the dimensionless inverse viscosity (Nf), Eötvös (Eo) and Froude numbers (Um), the latter being based on the centreline liquid velocity, using a Galerkin finite-element method. A boundary-fitted domain is used to examine the dependence of the steady bubble shape on a wide range of Nf, Eo and Um. Our analysis of the bubble nose and bottom curvatures shows that the intervals Eo=[20,30) and Nf=[60,80) are the limits below which surface tension and viscosity, respectively, have a strong influence on the bubble shape. In the interval Eo=(60,100], all bubble features studied are weakly dependent on surface tension. A linear stability analysis of the axisymmetric base states shows that there exist regions of (Nf,Eo,Um) space within which the bubble is unstable and assumes an asymmetric shape. To elucidate the mechanisms underlying the instability, an energy budget analysis is carried out which reveals that perturbation growth is driven by the bubble pressure for Eo≥100, and by the tangential interfacial stress for Eo<100. Examples of the asymmetric bubble shapes and their associated flow fields are also provided near the onset of instability for a wide range of Nf, Eo and Um.