Certain features of the Earth's magnetic field remain stationary for at least a sizeable fraction of the presumed convective time-scale of the outer core. These features may result from forced convection (laterally variable boundary buoyancy conditions) in the core as a result of core-mantle thermal coupling and the longer convective time-scale of the mantle. Although time-dependent motions in the core are certainly present, the observation of stationary features has motivated us to develop a method to find steady, finite-amplitude solutions to the non-linear convective equations. In this paper, we develop the iterative numerical method, and apply it to free convection (laterally uniform boundary buoyancy conditions), which is better understood than forced convection. Although we cannot yet be sure of the stability with respect to time perturbations of our solutions, our method has computational advantages over time-stepping that allow us to examine a large parameter space. We provide destabilization in our model by uniformly fluxing buoyancy across the bottom boundary with a concominant buoyancy sink at the top boundary. The Rayleigh number Ra is a measure of the destabilization. All calculations are at infinite Prandtl number Pr and magnetic Prandtl number q = 1. We first apply our method to a non-rotating, electrically insulating, spherical fluid shell, to demonstrate the method's viability. We then look for solutions in a rapidly rotating spherical fluid shell (high Taylor number; Ta = O(10 6)). We obtain axisymmetric polar modes, but not the non-axisymmetric equatorial modes (columns) predicted by linear theory, which require an azimuthal drift in spherical geometry. We next impose a poloidal magnetic field, and look for the magnetic analog of the polar modes. When the Lorentz force is comparable with the Coriolis force (El (Elsasser number) = O(10 0)), the modes fill the shell and most efficiently transport buoyancy. Thus, g z , the component of gravity parallel to the rotation axis, may play an important dynamical role for supercritical rotating convection, especially for an electrically conducting fluid in which the polar modes are more efficient than for an insulating fluid. In conjunction with the greater amplitude of convection near the inner spherical boundary than near the outer one, we find that the non-linear interaction between the advectively created toroidal magnetic field and the associated radial electrical current leads to a consistent bias towards equatorial upwelling flow. Until we determine the growth rates of instabilities, we cannot be sure of the stability of any particular converged solution, but we nevertheless believe that this non-linear result is a general feature of the ω-effect, with possible relevance to the poloidal circulation in the Earth's core.