The anharmonic phenomena in Zirconium Hydrides and Deuterides, including {epsilon}-ZrH2, {gamma}-ZrH, and {gamma}-ZrD, have been investigated from aspects of inelastic neutron scattering (INS) and lattice dynamics calculations within the framework of density functional theory (DFT). The observed multiple sharp peaks below harmonic multi-phonon bands in the experimental spectra of all three materials did not show up in the simulated INS spectra based on the harmonic approximation, indicating the existence of strong anharmonicity in those materials and the necessity of further explanations. We present a detailed study on the anharmonicity of zirconium hydrides/deuterides by exploring the 2D potential energy surface of hydrogen/deuterium atoms, and solving the corresponding 2D single-particle Schrodinger equation to get the eigenfrequencies. The obtained results well describe the experimental INS spectra and show harmonic behavior in the fundamental modes and strong anharmonicity at higher energies.