  A robust and general solver for the radial Schr\"odinger, Dirac, and
Kohn--Sham equations is presented. The formulation admits general potentials
and meshes: uniform, exponential, or other defined by nodal distribution and
derivative functions. For a given mesh type, convergence can be controlled
systematically by increasing the number of grid points. Radial integrations are
carried out using a combination of asymptotic forms, Runge-Kutta, and implicit
Adams methods. Eigenfunctions are determined by a combination of bisection and
perturbation methods for robustness and speed. An outward Poisson integration
is employed to increase accuracy in the core region, allowing absolute
accuracies of $10^{-8}$ Hartree to be attained for total energies of heavy
atoms such as uranium. Detailed convergence studies are presented and
computational parameters are provided to achieve accuracies commonly required
in practice. Comparisons to analytic and current-benchmark density-functional
results for atomic number $Z$ = 1--92 are presented, verifying and providing a
refinement to current benchmarks. An efficient, modular Fortran 95
implementation, \ttt{dftatom}, is provided as open source, including examples,
tests, and wrappers for interface to other languages; wherein particular
emphasis is placed on the independence (no global variables), reusability, and
generality of the individual routines.
