Spectral Density in the Gaussian Unitary Ensemble

In this post I perform numerical experiments on the spectral density in the Gaussian Unitary Ensemble (GUE).
This is the set of Hermitian $N \times N$ matrices $H$ with probability density $$\begin{equation*} P(H) dH = \mathcal{N} \prod_{ i < j} d \mathrm{Re} H_{ij} \prod_{ i < j} d \mathrm{Im} H_{ij}\ \prod_{ i} d H_{ii}\exp \left( -\frac{1}{2}\ \mathrm{tr}(H^2)\right) \end{equation*}$$ where $\mathcal{N}$ is the normalization constant. The spectral density $\rho (\lambda)$ is defined by $$\begin{equation*} \rho(x) dx = \sum_{k=1}^N P(x \le \lambda_k \le x + dx) \end{equation*}$$ with $\lambda_k$ the $N$ eigenvalues of the matrix $H$. In chapter 6 of his book, Mehta calculates that $\rho(x) = K_N(x,x)$ with $$\begin{equation*} K_N(x,y) = \frac{H_N(x) H_{N-1}(y) - H_N(y) H_{N-1}(x)}{x-y}e^{-\displaystyle\frac{x^2}{4} -\displaystyle\frac{y^2}{4}} \end{equation*}$$ Thus $$\begin{equation}\label{eq:20151114c} \rho(x) = \left( H_{N-1}(x) H_N'(x) - H_{N}(x) H_{N-1}'(x) \right) e^{-\displaystyle\frac{x^2}{2}} \end{equation}$$ Here $H_N(x)$ are the Hermite polynomials. I now check formula $\eqref{eq:20151114c}$ for $N = 1, 2, 3$ with a Monte Carlo simulation in Mathematica. I generate matrices from the GUE as follows: I draw matrices $A$ and $B$ whose elements $A_{ij}$ and $B_{ij}$ are independent normal random variables and take $$\begin{equation*} H = \frac{1}{2} \left( A + A^T \right) + \frac{i}{2} \left( B -B^T \right) \end{equation*}$$ Then I calculate the $N$ eigenvalues of $H$. This is the code that I use for this.

Clear[randomsamples] (* This code calculates the eigenvalues of NSamples matrices from GUE. The size of the matrices is NN times NN *) 
randomsamples[NN_, NSamples_] := Table[A = RandomVariate[NormalDistribution[], {NN, NN}]; 
B = RandomVariate[NormalDistribution[], {NN, NN}]; 
H = 1/2 (A + Transpose[A]) + I / 2 (B - Transpose[B]); 
Eigenvalues[H], {k, NSamples}] // Flatten

N = 1
In this case $$\begin{equation*} \rho(x) = \frac{1}{\sqrt{2 \pi }} \exp \left(-\frac{x^2}{2}\right) \end{equation*}$$ I produce a graph with the following commands  
NSamples = 10000; 
NN = 1; 
randomeigs = randomsamples[NN, NSamples]; 
fig1 = Histogram[randomeigs, {-3, 3, 0.1}, "PDF"]; 
fig2 = Plot[1/Sqrt[2 Pi] Exp[-x^2/2], {x, -3, 3}, PlotStyle -> Red]; 
Show[fig1, fig2]

N = 2
In this case $$\begin{equation*} \rho(x) = \frac{e^{-\frac{x^2}{2}} \left(x^2+1\right)}{2 \sqrt{2 \pi }} \end{equation*}$$ I produce a graph with the following commands
NSamples = 50000; NN = 2; 
randomeigs = randomsamples[NN, NSamples]; 
fig1 = Histogram[randomeigs, {-4, 4, 0.1}, "PDF"]; 
fig2 = Plot[ (E^(-(x^2/2)) (1 + x^2))/(2 Sqrt[2 \[Pi]]), {x, -4, 4}, PlotStyle -> Red]; 
Show[fig1, fig2]

N = 3
In this case $$\begin{equation*} \rho(x) = \frac{e^{-\frac{x^2}{2}} \left(x^4+3\right)}{6 \sqrt{2 \pi }} \end{equation*}$$ I produce a graph with the following commands
NSamples = 50000; NN = 3; 
randomeigs = randomsamples[NN, NSamples]; 
fig1 = Histogram[randomeigs, {-5, 5, 0.1}, "PDF"]; 
fig2 = Plot[ (E^(-(x^2/2)) (3 + x^4))/(6 Sqrt[2 \[Pi]]), {x, -5, 5}, PlotStyle -> Red]; Show[fig1, fig2]

The analytical expression $\eqref{eq:20151114c}$ thus agrees well with the Monte Carlo simulation.
Warning: In this post, I have not paid attention to overall constants in expressions, I only made sure the normalizations were consistent in the Mathematica code.