The possibility of bifurcation of the axisymmetric solution to non-axisymmetric buckling modes is examined by using a dedicated bifurcation check. The bifurcation check was found to yield such buckling modes correctly. The evolution of such buckling modes is simulated by a separate non-axisymmetric model assigned imperfections with axisymmetric and nonaxisymmetric components. The domain analyzed is one characteristic wavelength long (2λ_C). Initially, compression activates mainly axisymmetric deformation. In the neighborhood of the bifurcation point, non-axisymmetric deformation starts to develop, eventually leading to a limit load instability. Experimental responses were simulated with accuracy by assigning appropriate values to the two imperfection amplitudes. Prediction of the limit strains for the whole range of diameter-to-thickness ratios (D/t) considered in the experiments was achieved by making the amplitude of the non-axisymmetric imperfection proportional to (D/t)²/m³ (m is the circumferential wavenumber). Matching all aspects of the experiments required inclusion of the anisotropy measured in the tubes tested through Hill's yield criterion in all models.