AbstractThe nodal integral method is a relatively new numerical technique that has been used recently to solve both static and dynamic multidimensional problems in heat transfer, fluid flow and neutron transport. The method offers significant advantages in terms of stability, accuracy and efficiency over conventional finite elements when the problem can be adequately modelled in Cartesian co‐ordinates. This method was used to investigate bifurcation phenomena in the Bénard problem for aspect ratios in the range of one to nine. Automatic search techniques were used with a static version to find the first four critical Rayleigh values for a square cavity, to map the first two critical Rayleigh values as a function of aspect ratio, and to examine the solution types. Accuracy enhancement was obtained by factorization and extrapolation. Critical values, obtained by interpolation, were verified dynamically. Aspect ratio crossover and transition values were found for the first two critical Rayleigh numbers, with an accuracy of the order of ±3 per cent. The precision achieved in the results for Ra* and Ra** as a function of β is usually within 0.1%–0.2% except at high β (i.e. near β=9.0) and at large critical values of Ra (i.e. the first few values of Ra** near β=1). Specific results at β=1.0 are Ra*=2584±0.5, Ra**=6807, Ra3* = 19 734 and Ra4*=22 586.