Thin-walled shells of revolution under circumferentially uniform pre-buckling stress states are important fundamental systems, often serving as reference 'base cases' to which the behaviour of more complex unsymmetrical systems can be related. However, the same simplicity that often permits closed-form algebraic expressions for the critical buckling load is also often responsible for a lack of localisation and significant ambiguity in the critical buckling mode. The computation of the linear or nonlinear buckling loads requires the systematic trial of many potential buckling modes to identify the one which minimises the necessary strain energy. In times when researchers used custom-written tools usually based on circumferential Fourier series expansions, this operation was relatively straightforward. However, today's analysts using 'general' commercial 3D finite element packages must apply careful safeguards to correctly identify the correct buckling load and associated mode in axisymmetric shell systems. This paper presents a detailed computational strategy to accurately and efficiently investigate the correct buckling load and mode number of axisymmetric shell systems using the technique of a 'panel analysis'. This is implemented in the ABAQUS finite element solver controlled by the SIMULIA TM Isight automation software and the Python object-oriented programming language. The methodology is illustrated on three classical benchmark problems from the scientific literature on the buckling of cylindrical shells under meridional compression, with special attention given to meshing considerations.