Background
Let $M$ be a $N\times N$ random symmetric matrix, which is drawn from a Gaussian distribution of mean $\mu=0$ and variance $\sigma^2=1/N$.
The first level correlation is defined by: $$\begin{equation} \rho(\lambda)=\left \langle \text{Tr}\delta(\lambda-M)\right \rangle\quad\quad\quad\quad \tag{1} \label{eqn1} \end{equation}$$ Which is nothing more than the distribution of eigenvalues: following my routine below I can compute this quantity. The brackets represents the average over the gaussian distribution.
RandomMatrix[n_] :=
RandomMatrix[n] =
RandomVariate[NormalDistribution[0, 1/Sqrt[n]], {n, n}];
n := 1500
M := 1/Sqrt[2]*(RandomMatrix[n] + Transpose@RandomMatrix[n])
\[Rho] = Eigenvalues[M];
Histogram[\[Rho]]
In the large limit $N\gg 1$ equation \eqref{eqn1} will describe a semi-circular law.
Question
The two-level correlation function is defined by: $$\begin{equation} \rho^{(2)}(\lambda, \mu)=\left\langle\frac{1}{N} \operatorname{Tr} \delta(\lambda-M) \frac{1}{N} \operatorname{Tr} \delta(\mu-M)\right\rangle\quad\quad\quad\quad \tag{2} \label{eqn2} \end{equation}$$ I interpret this quantity is as how correlated two given eigenvalues will be. Following the previous routine, my $\rho$ will be a list of the $N$ eigenvalues of $M$. How can I numerically compute equation \eqref{eqn2}? Isn't it impossible to tell how correlated two elements are in one same list? If I generate several lists, it has to be for different matrices, which does not work anymore.
Given one matrix $M$, how can I compute the correlation between two eigenvalues?
Any help or remark is appreciated!

RandomMatrix[n]to stay the same when I call it twice:RandomMatrix[n]+Transpose@RandomMatrix[n]. The second reason was to emphasize the fact that when I speak of computing the correlation function, it is for a single matrix $M$ and not a collection of random matrices. For example when I computed the eigenvalues it was for one matrix, I didn't average it over many matrices. $\endgroup$RandomMatrixdon't have the same distribution as the off-diagonal element (Variance[Diagonal[M]] / Variance[Flatten[M]]is about 2). The off-diagonals are the sum of 2 independent RVs, while the diagonal is the sum of the same RV. $\endgroup$