We develop systematic formulations for calculating the magnetization of spins diffusing in a bounded region in the presence of the surface relaxation and magnetic field inhomogeneity and compute explicitly the relaxation exponent for the Carr-Purcell-Meiboom-Gill spin echoes. The results depend on the echo number n, and three dimensionless parameters: L-rho/L-S, (D) over tilde (0) =(L-D/L-S)(2), the dimensionless diffusion constant, and <()over tilde>=(LDLS)-L-2/L-G(3)=Delta omega tau, the dimensionless gyromagnetic ratio, where the restriction is characterized by a size L-S, the magnetic field inhomogeneity by a dephasing length, L-G, the diffusion length during half-echo time by L-D, and a length L-rho characterizes the surface relaxation. Here Delta omega is the line broadening and 2 tau is the echo period. Depending on the length scales, three main regimes of decay have been identified: short-time, localization, and motionally averaging regimes (MAv). The short-time and the MAv regimes are described well by the cumulant expansion in terms of powers of the "small" parameter <()over tilde>. We give simple formulas for decay rates in these two asymptotic regimes. We show that the Gaussian phase approximation (GPA), i.e., the exponent up to the second order in <()over tilde>(2) in terms of a full eigenmode expansion interpolates well between these two regimes. In the localization regime, the decay exponent depends on a fractional power, <()over tilde>(2/3), denoting a breakdown of the GPA and a breakdown of the cumulant expansion in terms of <()over tilde>.