The nearest-neighbor level spacing distribution is numerically investigated by directly diagonalizing disordered Anderson Hamiltonians for systems of sizes up to 100×100×100 lattice sites. The scaling behavior of the level statistics is examined for large spacings near the delocalization-localization transition and the correlation length exponent is found. By using high-precision calculations we conjecture a new interpolation of the critical cumulative probability, which has size-independent asymptotic form ln I(s) ∝ −s α with α = 1.0 ± 0.1. The statistical fluctuations in energy spectra of disordered quantum systems attract at present much attention [1,2,3,4,5]. It is known that by increasing the fluctuations of a random potential the one-electron states undergo a localization transition, which is the origin of the Anderson metal-insulator transition (MIT) [6]. The influence of the disorder on the wave functions is reflected by the mutual correlations between the corresponding energy levels, so that the statistics of energy levels is sensitive to the MIT. In the metallic limit the statistics of energy spectra can be described by the random-matrix theory (RMT) developed by Wigner and Dyson [7,8]. This was shown by solving the zero-mode nonlinear σ-model using the supersymmetric formalism [9]. Later, perturbative corrections to the two-level correlation function obtained in the RMT were evaluated in the diffusive regime by the impurity diagram technique [10]. In the insulating regime, when the degree of disorder W is much larger than the critical value W c , the energy levels of the strongly localized eigenstates fluctuate as independent random variables.An important quantity for analyzing the spectral fluctuations is the nearest-neighbor level spacing distribution P (s). It contains information about all of the n−level correlations. In the metallic regime P (s) is very close to the Wigner surmise P W (s) = π s/2 exp −π s 2 /4 [11] (s is measured in units of the mean level spacing ∆). In the localized regime the spacings are distributed according to the Poisson law, P P (s) = exp(−s), because the levels are completely uncorrelated. The study of the crossover of P (s) between the Wigner and the Poissonian limits which accompanies the disorder-induced MIT in threedimensional system (3D) was started in Ref. [12] and became the subject of several subsequent investigations [1,2,5,13,14].It was suggested earlier [1] that P (s) exhibits critical behavior and should be size-independent at the MIT. Investigating the finite-size scaling properties of P (s) provides not only an alternative method for locating the transition [2], but allows also to determine the critical behavior of the correlation length [14]. A technical advantage of the method is that one needs to compute only energy spectra and not eigenfunctions and/or the conductivity. On the other hand, a large number of realizations of the random potential has to be considered. In comparison with the well-established transfer-matrix method [15] by which one approache...