Morphing shells are nonlinear structures that have the ability to change shape and adopt multiple stable states. By exploiting the concept of morphing, designers may devise adaptable structures, capable of accommodating a wide range of service conditions, minimising design complexity and cost. At present, models predicting shell multistability are often characterised by a compromise between computational efficiency and result accuracy. This paper addresses the main challenges of describing the multistable behaviour of thin composite shells, such as bifurcation points and snap-through loads, through the development of an accurate and computationally efficient energy-based method. The membrane and the bending components of the total strain energy are decoupled by using the semi-inverse formulation of the constitutive equations. Transverse displacements are approximated by using Legendre polynomials and the membrane problem is solved in isolation by combining compatibility conditions and equilibrium equations. This approach provides the strain energy as a function of curvature only, which is of particular interest, as this decoupled representation facilitates efficient solution. The minima of the energy with respect to the curvature components give the multiple stable configurations of the shell. The accurate evaluation of the membrane energy is a key step in order to correctly capture the multiple configurations of the structure. Here, the membrane problem is solved by adopting the Differential Quadrature Method (DQM), which provides accurate results at a relatively small computational cost. The model is benchmarked against three exemplar case studies taken from the literature.