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.