The stability of two-dimensional (2D) natural convection flows with respect to both two- and three-dimensional perturbations is investigated numerically. Several methods (Arnoldi’s method, preconditioned Newton’s iteration and preconditioned continuation method) are put together for this purpose and applied to natural convection in a differentially heated square cavity with conducting horizontal walls for a large range of Prandtl numbers. These methods are first validated by comparison with results reported in the literature for several Prandtl numbers for both two-dimensional and three-dimensional perturbations. They are then used to extend the stability of the 2D base flows to two-dimensional perturbations for other Prandtl numbers, and to investigate in detail the stability of these two-dimensional base solutions with respect to three-dimensional perturbations for a large range of Prandtl number. For Prandtl number equal to 1 the base solutions are more unstable to three-dimensional perturbations and the most unstable three-dimensional mode, which is oscillatory, is connected to the third most unstable two-dimensional mode; for larger Prandtl numbers (3, 7, and 20), the 2D base flows are found more stable to three-dimensional perturbations than to the 2D ones and the most unstable modes are thus two dimensional; for smaller Prandtl numbers (0.71, 0.1, and 0.015) the base solutions are again more unstable to three-dimensional perturbations. These three-dimensional modes are however no longer connected to the most unstable two-dimensional oscillatory modes and correspond to stationary perturbations. Furthermore for Pr=0.015 and 0.1 the critical Rayleigh number is approximately 40 times smaller than that found for two-dimensional perturbations, indicating a dramatic change in the nature of the instability which seems to be linked to Taylor–Görtler-type instability.