This paper proposes an effective diffusion equation method to analyze nuclear magnetic resonance (NMR) relaxation. NMR relaxation is a spin system recovery process, where the evolution of the spin system is affected by the random field due to Hamiltonians, such as dipolar couplings. The evolution of magnetization can be treated as a random walk in phase space described either by a normal or fractional phase diffusion equation. Based on these phase diffusion equations, the NMR relaxation rates and equations can be obtained, exemplified in the analysis of relaxations affected by an arbitrary random field, and by dipolar coupling for both like and unlike spins. The obtained theoretical results are consistent with the reported results in the literature. Additionally, the anomalous relaxation expression obtained from the Mittag-Leffler function based time correlation function can successfully fit the previously reported 13C T1 NMR experimental data of polyisobutylene (PIB) in the blend of PIB and head-to-head poly(propylene) (hhPP). Furthermore, the proposed phase diffusion approach provides an intuitive way to interpret NMR relaxation, particularly for the fractional NMR relaxation, which is still a challenge to explain by the available theoretical methods. The paper provides additional insights into NMR and magnetic resonance imaging (MRI) relaxation experiments.