  We report on a numerical study of the density matrix functional introduced by
Lieb, Solovej and Yngvason for the investigation of heavy atoms in high
magnetic fields. This functional describes {\em exactly} the quantum mechanical
ground state of atoms and ions in the limit when the nuclear charge $Z$ and the
electron number $N$ tend to infinity with $N/Z$ fixed, and the magnetic field
$B$ tends to infinity in such a way that $B/Z^{4/3}\to\infty$. We have
calculated electronic density profiles and ground state energies for values of
the parameters that prevail on neutron star surfaces and compared them with
results obtained by other methods. For iron at $B=10^{12}$ G the ground state
energy differs by less than 2 \% from the Hartree-Fock value. We have also
studied the maximal negative ionization of heavy atoms in this model at various
field strengths. In contrast to Thomas-Fermi type theories atoms can bind
excess negative charge in the density matrix model. For iron at $B=10^{12}$ G
the maximal excess charge in this model corresponds to about one electron.
