A new approximate method is presented for calculating the density of states and one-electron Green's function in the low-energy tail of a high-density impurity band. Such states occur in regions of large attractive potential produced by unusually high (random) concentrations of attractive centers and/or unusually low concentrations of repulsive centers. These well-separated deep wells each have one bound state of lowest energy. The distribution of these lowest levels vields the low-energy tail. For any one well, a bound-state energy $E({\mathbf{x}}_{0})$ is estimated variationally using $\ensuremath{\psi}(\mathbf{x})=f(\mathbf{x}\ensuremath{-}{\mathbf{x}}_{0})$, where $f(\mathbf{x})$ has any fixed form. The best energy for any well is obtained by minimizing $E({\mathbf{x}}_{0})$ with respect to ${\mathbf{x}}_{0}$ in the vicinity of that well. The number of wells that contribute to the density of states at energy $E$ is given by the number of local minima in $E({\mathbf{x}}_{0})$ that occur at the level $E$. In the high-density limit, Gaussian statistics are adequate for treating the potential fluctuations and the expectation value for the density of states is easily calculated. At low energies the best choice for the function $f$ is that which maximizes the expected density of states. Application to an exactly soluble one-dimensional model, a Gaussian white noise potential, yields the correct asymptotic form $\ensuremath{\rho}(E)=\mathrm{const}|E|\mathrm{exp}[\ensuremath{-}(\frac{4}{3}){|2E|}^{\frac{3}{2}}]$. In three dimensions, the density of states in the Gaussian approximation is found to have the form $\ensuremath{\rho}(E)=[\frac{A(E)}{{\ensuremath{\xi}}^{2}}]\mathrm{exp}[\frac{\ensuremath{-}B(E)}{(2\ensuremath{\xi})}]$, where $\ensuremath{\xi}$ is proportional to the concentration of impurities and to the square of the strength of the impurity potential. For screened, charged impurities, $B(E)={E}_{{Q}^{2}}b(\ensuremath{\nu})$, $A(E)={Q}^{3}{E}_{{Q}^{3}}a(\ensuremath{\nu})$, where ${E}_{Q}=\frac{{\ensuremath{\hbar}}^{2}{Q}^{2}}{(2{m}^{*})}$, ($\frac{1}{Q}$) is the screening radius, and $\ensuremath{\nu}=\frac{({E}_{0}\ensuremath{-}E)}{{E}_{Q}}$ is the energy below the mean potential ${E}_{0}$ in units of ${E}_{Q}$. Computer-calculated curves are provided for the universal dimensionless functions $a(\ensuremath{\nu})$ and $b(\ensuremath{\nu})$. The exponent $b(\ensuremath{\nu})$ behaves like ${\ensuremath{\nu}}^{\frac{1}{2}}$ when $\ensuremath{\nu}$ is small (strong screening) and behaves like ${\ensuremath{\nu}}^{2}$ when $\ensuremath{\nu}$ is large (weak screening).