Needs["Histograms`"]

(*Mathematica example for random matrix exercise*)(*Generating M members of the GOE ensemble of NxN matrices,and calculating the difference between the eigenvalues in the middle of the range*)

(*GENERATING RANDOM GOE MATRICES*)

(*Generating M members of the GOE ensemble of NxN matrices,and calculating the difference between the eigenvalues in the middle of the range*)

(*Increase M,N as appropriate*)

M=20;
Size=4;

(*Generate M random SizexSize arrays of Gaussian (normal) random numbers;symmetrize them;find their sorted eigenvalues,keep a list of the difference between the two middle ones*)diffs=Table[0,{i,1,M}];
M11=Table[0,{i,1,M}];
For[i=1,i≤M,i++,{Mat=Table[Random[NormalDistribution[]],{i,1,Size},{j,1,Size}];
MatSym=Mat+Transpose[Mat];
eigs=Sort[Eigenvalues[MatSym]];
diffs[[i]]=eigs[[Size/2+1]]-eigs[[Size/2]]};
  M11[[i]]=MatSym[[1]][[1]];]

(*You can also do this with functional programming:it's considered more elegant in Mathematica world*)Symmetrize[m_]:=m+Transpose[m];
EigSorted[m_]:=Sort[Eigenvalues[m]]
EigDiffs[eigs_]:=eigs[[Size/2+1]]-eigs[[Size/2]]
diffs=Map[EigDiffs[EigSorted[Symmetrize[#]]]&,Table[Random[NormalDistribution[]],{m,1,M},{i,1,Size},{j,1,Size}]];
M11 = Map[Symmetrize[#][[1,1]]&,Table[Random[NormalDistribution[]],{m,1,M},{i,1,Size},{j,1,Size}]];

(*Divide out by mean value of the splittings*)
diffAve=Mean[diffs]

(*Histogram generates a histogram of diffs.HistogramScale normalizes the integral to one.*)Histogram[diffs/diffAve,HistogramScale->1]

(*Demonstrated with histogram for diagonal element First plot histogram for M11;store as "histGauss" for later plotting with theory*)
histGauss=Histogram[M11,HistogramScale->1]

(*Now,plot expected Gaussian fit*)sigma=2;
theoryGauss=Plot[(1/Sqrt[2 Pi sigma^2]) Exp[-x^2/(2 sigma^2)],{x,-5,5}]

(*and show the comparison*)

Show[histGauss,theoryGauss]

(*GENERATING RANDOM+-1 MATRICES*)

(*Symmetric matrix with integer values+-1 with 50/50 probability*)

(*Start with random matrix of+-1,no symmetry*)

MatPM=Table[Random[Integer,{0,1}]*2-1,{i,1,Size},{j,1,Size}];

MatPMSym=Table[If[i≥j,MatPM[[i]][[j]],MatPM[[j]][[i]]],{i,1,Size},{j,1,Size}];

(*Take a look at it*)
MatrixForm[MatPMSym]