  We introduce \texttt{featom}, an open source code that implements a
high-order finite element solver for the radial Schr\"odinger, Dirac, and
Kohn-Sham equations. The formulation accommodates various mesh types, such as
uniform or exponential, and the convergence can be systematically controlled by
increasing the number and/or polynomial order of the finite element basis
functions. The Dirac equation is solved using a squared Hamiltonian approach to
eliminate spurious states. To address the slow convergence of the $\kappa=\pm1$
states due to divergent derivatives at the origin, we incorporate known
asymptotic forms into the solutions. We achieve a high level of accuracy
($10^{-8}$ Hartree) for total energies and eigenvalues of heavy atoms such as
uranium in both Schr\"odinger and Dirac Kohn-Sham solutions. We provide
detailed convergence studies and computational parameters required to attain
commonly required accuracies. Finally, we compare our results with known
analytic results as well as the results of other methods. In particular, we
calculate benchmark results for atomic numbers ($Z$) from 1 to 92, verifying
current benchmarks. We demonstrate significant speedup compared to the
state-of-the-art shooting solver \texttt{dftatom}. An efficient, modular
Fortran 2008 implementation, is provided under an open source, permissive
license, including examples and tests, wherein particular emphasis is placed on
the independence (no global variables), reusability, and generality of the
individual routines.
