I tried to use Mathematica to study random matrix and expect to get Distribution of level spacings $p_2(s)=\frac{32}{\pi^2}s^2 \exp(-\frac{4}{\pi}s^2)$ for large Gaussian unitary matrix. However, the PDF and CDF histograms did not fit very well. It is NOT due to random fluctuation. Is there something wrong with me?
n = 4000;
a = RandomVariate[GaussianUnitaryMatrixDistribution[n]];
s = Eigenvalues[a] // Sort;
s = Subsequences[s, {2}];
t = (#[[2]] - #[[1]]) & /@ s;
t = t/Mean[t];
p2[x_] := 32/\[Pi]^2 x^2 E^(-4/\[Pi] x^2);
Show[{Histogram[t, {0, 3, 0.01}, "PDF"],
Plot[p2[x], {x, 0, 3}]}]
Show[{Histogram[t, {0, 3, 0.01}, "CDF"],
Plot[Evaluate[Integrate[p2[x], {x, 0, x}]], {x, 0, 3}]}]
