A flexible Galerkin method based on proper orthogonal decomposition (POD) is described to construct the bifurcation diagram, as the Rayleigh number R is varied, in the Rayleigh–Bénard convection in a rectangular box for large Prandtl number. The bifurcation diagram is approximated using the POD modes resulting from unconverged snapshots for just one specific value of R, calculated in either Newton iterations or time-dependent runs converging to steady states. Moreover, the selection of the specific value of R is quite flexible. In addition, a horizontal reflection symmetry is taken into account to construct a symmetry-preserving Galerkin system. The resulting un-symmetric and symmetric low-dimensional systems are combined with a basic continuation method, which provide the bifurcation diagram at a quite low computational cost.