Earth’s magnetic field is generated in its fluid metallic core through motional induction in a process termed the geodynamo. Fluid flow is heavily influenced by a combination of rapid rotation (Coriolis forces), Lorentz forces (from the interaction of electrical currents and magnetic fields) and buoyancy; it is believed that the inertial force and the viscous force are negligible. Direct approaches to this regime are far beyond the reach of modern high-performance computing power, hence an alternative ‘reduced’ approach may be beneficial. Taylor (Proc. R. Soc. Lond. A, vol. 274 (1357), 1963, pp. 274–283) studied an inertia-free and viscosity-free model as an asymptotic limit of such a rapidly rotating system. In this theoretical limit, the velocity and the magnetic field organize themselves in a special manner, such that the Lorentz torque acting on every geostrophic cylinder is zero, a property referred to as Taylor’s constraint. Moreover, the flow is instantaneously and uniquely determined by the buoyancy and the magnetic field. In order to find solutions to this mathematical system of equations in a full sphere, we use methods of optimal control to ensure that the required conditions on the geostrophic cylinders are satisfied at all times, through a conventional time-stepping procedure that implements the constraints at the end of each time step. A derivative-based approach is used to discover the correct geostrophic flow required so that the constraints are always satisfied. We report a new quantity, termed the Taylicity, that measures the adherence to Taylor’s constraint by analysing squared Lorentz torques, normalized by the squared energy in the magnetic field, over the entire core. Neglecting buoyancy, we solve the equations in a full sphere and seek axisymmetric solutions to the equations; we invoke $\unicode[STIX]{x1D6FC}$- and $\unicode[STIX]{x1D714}$-effects in order to sidestep Cowling’s anti-dynamo theorem so that the dynamo system possesses non-trivial solutions. Our methodology draws heavily on the use of fully spectral expansions for all divergenceless vector fields. We employ five special Galerkin polynomial bases in radius such that the boundary conditions are honoured by each member of the basis set, whilst satisfying an orthogonality relation defined in terms of energies. We demonstrate via numerous examples that there are stable solutions to the equations that possess a rapidly decreasing spectrum and are thus well-converged. Classic distributions for the $\unicode[STIX]{x1D6FC}$- and $\unicode[STIX]{x1D714}$-effects are invoked, as well as new distributions. One such new $\unicode[STIX]{x1D6FC}$-effect model possesses oscillatory solutions for the magnetic field, rarely before seen. By comparing our Taylor state model with one that allows torsional oscillations to develop and decay, we show the equilibrium state of both configurations to be coincident. In all our models, the geostrophic flow dominates the ageostrophic flow. Our work corroborates some results previously reported byWu & Roberts (Geophys. Astrophys. Fluid Dyn., vol. 109 (1), 2015, pp. 84–110), as well as presenting new results; it sets the stage for a three-dimensional implementation where the system is driven by, for example, thermal convection.