The Grad–Shafranov equation (GSE) for axisymmetric MHD equilibria is a nonlinear, scalar PDE which in principle can have zero, one or more non-trivial solutions. The conditions for the existence of multiple solutions has been little explored in the literature so far. We develop a simple analytic model to calculate multiple solutions in the large aspect ratio limit. We compare the results to the recently developed deflated continuation method to find multiple solutions in a realistic geometry and right-hand side of the GSE using the finite element method. The analytic model is surprisingly accurate in calculating multiple solutions of the GSE for given boundary conditions and the two methods agree well in limiting cases. We examine the effect of plasma shaping and aspect ratio on the multiple solutions and show that shaping generally does not alter the number of solutions. We discuss implications for predictive modelling, equilibrium reconstruction, plasma stability and disruptions.