A number of computational schemes for the calculation of N-15 and P-31 NMR chemical shifts and shielding constants in a series of azoles, phospholes, and phosphazoles was examined. A very good correlation between calculated at the CCSD(T) level and experimental N-15 and P-31 NMR chemical shifts was observed. It was found that basically solvent, vibrational, and relativistic corrections are of the same order of magnitude and alternate in sign, being, on average, of about 2-3 ppm in absolute value but, being much larger (up to 14 ppm) in the case of solvent molecules explicitly introduced into computational space. At the DFT level, the performance of nine exchange correlation functionals including six conventional gradient functionals and three hybrid functionals was studied. The most accurate results were reached with the OLYP and Keal-Tozer's family of functionals, KT1, KT2, and KT3, while the most popular B3LYP and PBE0 functionals showed the most unreliable results. On the basis of these data, we highly recommend OLYP and KT2 functionals for the computation of N-15 and P-31 NMR chemical shifts at the DFT level in the diverse series of nitrogen- and phosphorus-containing heterocycles. Benchmark calculations of N-15 and P-31 NMR chemical shifts in a series of larger nitrogen and phosphorus-containing heterocycles were performed at the DFT level in comparison with experiment and revealed the OLYP functional in combination with the aug-pcS-3/aug-pcS-2 locally dense basis set scheme as the most effective computational scheme.