Using the CASSCF/CASPT2 approach, along with several DFT methods (PBE0, B3LYP, BP86, OLYP), we have investigated the bonding of CO, NO, and O-2 molecules to two model heme systems: an iron(II) porphyrin with and without an axial imidazole ligand. The experimentally available binding energies are best reproduced by the CASPT2 method and with the OLYP functional. The other functionals considered perform much worse, either severly overbinding (BP86) or underbinding (B3LYP, PBE0). Significant discrepancies between the different density functionals are observed, not only for the energetics but sometimes also for structure predictions. This confirms our viewpoint that a balanced treatment of the electronic exchange and correlation is vital to describe the weak metal-ligand bond between heme and CO, NO, or O-2. The binding energies Delta E-b were split into two contributions: the so-called spin-pairing energy Delta E-sp and the "inherent" binding energy Delta E-b0, and both contributions were analyzed in terms of method and basis set effects. We have also investigated the spin density distributions resulting from the bonding of the NO molecule (a noninnocent ligand) to heme. Our analysis at the DFT and CASSCF level shows that, while various density functionals predict qualitatively very different spin distributions, the CASSCF spin populations most closely correspond to the results obtained with the pure BP86 or OLYP rather than with the hybrid functionals.