Multiple steady states are investigated for natural convection of fluids in a square enclosure with non-isothermally hot bottom wall, isothermally cold side walls, and thermally insulated top wall. A robust computation scheme involving steady-state governing equations has been developed to compute the steady states as a function of Rayleigh number (10 2 ≤ R a ≤ 10 6) for two different Prandtl numbers (Pr = 0.026 and 0.1). Penalty Galerkin finite element method with Newton-Raphson solver is employed for the solution of the governing equations, while the solution branches are initiated by varying initial guess to the Newton-Raphson solver. In this context, a dual-perturbation scheme involving perturbations of the boundary conditions and various process parameters has been designed leading to the rich spectrum of the symmetric and asymmetric solution branches for the current symmetric problem. It is found that multiple steady states occur beyond a critical value of Ra, which depends on the magnitude of Pr. In addition to the basic solution branch (corresponding to the solutions obtained via uniform initial guesses), nineteen additional solution branches (six symmetric and thirteen asymmetric) are obtained for Pr = 0.026, while four additional solution branches (two symmetric and two asymmetric) are obtained for Pr = 0.1. The solution branches are associated with a wide spectrum of flow structures (24 distinct types for Pr = 0.026 excluding the reflection symmetric mirror images of the asymmetric solutions), which are reported for the first time. The flow structures lead to various heating scenarios within the enclosure resulting in a significant variation of heat transfer rates (more than 50%). The current results are important for the practical applications. The spectrum of the possible scenarios revealed in this work can be pivotal to design the optimal processes based on the process requirement (targeted heating or enhanced heating rates). © 2021 Author(s).