Density functional theory (DFT) calculations on iron(II) (FeP) and cobalt(II) (CoP) porphyrins with nitric oxide (NO) and carbon monoxide (CO) as axial ligands have been performed using BLYP, B3LYP, OLYP, B3PW91, and M06L functionals with double- and triple-ξ basis sets. Optimized geometries and binding energies were found to depend very strongly on the functional and basis set used. MP2 fails completely to describe geometries and binding energies of FeP–CO, CoP–CO, and CoP–NO complexes, but performs relatively well for the FeP–NO complex. The calculated binding energies range from −4.5 to −19.6 kcal mol−1 for CoP–CO, +4.0 to −29.1 kcal mol−1 for FeP–CO, +7.6 to −41.1 kcal mol−1 for FeP–NO, and −2.2 to −36.7 kcal mol−1 for CoP–NO. B3PW91/6–311 + G(d,p) and OLYP/6–311+G(d,p) perform better than the other functional/basis set combinations in reproducing the available experimental data. A complete active space self-consistent field/complete active space with second-order perturbation theory (CASSCF/CASPT2) analysis of the CO and NO interactions with FeP and CoP can be interpreted either in terms of σ-donation/π-back-donation mechanism (FeP–CO and CoP–CO) or σ bond formation between the orbital of metal and the π* orbital of NO.