This work reports results of numerical simulations of viscous incompressible flow past a sphere. The primary objective is to identify transitions that occur with increasing Reynolds number, as well as their underlying physical mechanisms. The numerical method used is a mixed spectral element/Fourier spectral method developed for applications involving both Cartesian and cylindrical coordinates. In cylindrical coordinates, a formulation, based on special Jacobi-type polynomials, is used close to the axis of symmetry for the efficient treatment of the ‘pole’ problem. Spectral convergence and accuracy of the numerical formulation are verified. Many of the computations reported here were performed on parallel computers. It was found that the first transition of the flow past a sphere is a linear one and leads to a three-dimensional steady flow field with planar symmetry, i.e. it is of the ‘exchange of stability’ type, consistent with experimental observations on falling spheres and linear stability analysis results. The second transition leads to a single-frequency periodic flow with vortex shedding, which maintains the planar symmetry observed at lower Reynolds number. As the Reynolds number increases further, the planar symmetry is lost and the flow reaches a chaotic state. Small scales are first introduced in the flow by Kelvin–Helmholtz instability of the separating cylindrical shear layer; this shear layer instability is present even after the wake is rendered turbulent.