function z=levels(n,im)

%Plot exact semicircle formula
%The formula (5.2.16) in Mehta is the sum of phi_j^2
% but Christoffel-Darboux seems smarter

x=[-1:.001:1]*sqrt(2*n)*1.3;
pold = 0*x;                      % -1st Hermite polynomial
p= 1+0*x;                        %  0th Hermite polynomial
k=p;

for j=1:n;                       % Three term recurrence
   pnew = (sqrt(2)*x.*p-sqrt(j-1)*pold)/sqrt(j);
   pold = p; p=pnew;
end
pnew = (sqrt(2)*x.*p-sqrt(n)*pold)/sqrt(n+1);

k = n*p.^2-sqrt(n*(n+1))*pnew.*pold;% Use p.420 of Mehta

if (~im),                         % Extra term on p. 128
  k=k-sqrt(n/2)*pold.*(2*x.^2+1)/sqrt(3);
end

k=k.*exp(-x.^2)/sqrt(pi);

if ((~im)&rem(n,2))                    % alpha term on p.128
  k=k+ pold*gamma((n+1)/2)*2^(n/2-1)./(sqrt(pi*gamma(n))*exp(x.^2/2));
end

plot(x/sqrt(2*n),k*pi/sqrt(2*n),'b'); % Rescale so that "semicircle"
                                      % is on [-1,1] and area is pi/2

