<p>We develop and analyze a robust Arrow–Hurwicz (AH) iterative method for the numerical solution of steady Boussinesq flows with nonhomogeneous partitioned Dirichlet boundary conditions. Although a direct AH formulation may be applied to the momentum equation alone, we demonstrate that incorporating an AH-type update for the temperature equation is crucial for stability and convergence in buoyancy-driven systems, particularly at high Rayleigh numbers. The resulting Improved Arrow–Hurwicz (IAH) scheme avoids solving saddle-point systems at each iteration and yields a fully decoupled algorithm with low computational cost per step. We establish existence, uniqueness, uniform boundedness, and convergence under standard small-data assumptions, and provide corresponding error estimates for the finite element discretization. Extensive two- and three-dimensional numerical experiments verify the theoretical findings, demonstrate significant acceleration over the alternative AH scheme and the Penalty–Picard iteration, and confirm robust convergence in high–Rayleigh number regimes. The proposed method offers a scalable and efficient solver for steady natural convection and provides a promising alternative to continuation-based approaches traditionally used for high–Rayleigh flows.</p>