The coupled diffusion equations recently proposed by Tokuyama for concentrated hard-sphere suspensions are numerically solved, starting from nonequilibrium initial configurations. The most important feature of those equations is that the self-diffusion coefficient Ds(Φ) becomes zero at the glass transition volume fraction φg as Ds(Φ) approx. D0|1 - Φ(x,t)/φg|γ with γ = 2 where Φ(x,t) is the local volume fraction of colloids, D0 the single-particle diffusion constant, and φg = (4/3)3/(7 ln 3 - 8 ln 2 + 2). This dynamic anomaly results from the many-body correlations due to the long-range hydrodynamic interactions. Then, it is shown how small initial disturbances can be enhanced by this anomaly near φg, leading to long-lived, spatial heterogeneities. Those heterogeneities are responsible for the slow relaxation of nonequilibrium density fluctuations. In fact, the self-intermediate scattering function is shown to obey a two-step relaxation around the β-relaxation time tβ approx. |1 - φ/φg|-1, and also to be well approximated by the Kohlrausch-Williams-Watts function with an exponent β around the α-relaxation time tα approx. |1 - φ/φg|-η, where η = γ/β, and φ is the particle volume fraction. Thus, the nonexponential α relaxation is shown to be explained by the existence of long-lived, spatial heterogeneities.