We introduce a new computational method to study porphyrin-like transition metal complexes, bridging density functional theory and exact many-body techniques, such as the density matrix renormalization group (DMRG). We first derive a multi-orbital Anderson impurity Hamiltonian starting from first principles considerations that qualitatively reproduce GGA+U results when ignoring inter-orbital Coulomb repulsion $U$ and Hund exchange $J$. An exact canonical transformation is used to reduce the dimensionality of the problem and make it amenable to DMRG calculations, including all many-body terms (both intra, and inter-orbital), which are treated in a numerically exact way. We apply this technique to FeN$_4$ centers in graphene and show that the inclusion of these terms has dramatic effects: as the iron orbitals become single occupied due to the Coulomb repulsion, the inter-orbital interaction further reduces the occupation yielding a non-monotonic behavior of the magnetic moment as a function of the interactions, with maximum polarization only in a small window at intermediate values of the parameters. Furthermore, $U$ changes the relative position of the peaks in the density of states, particularly on the iron $d_{z^2}$ orbital, which is expected to greatly affect the binding of ligands.