Function acts bizarrely

Let's define X as :

and objects releated to it:

I want to calculation value of function following, and plot it on the same graph with eigenvalues of Y.

My work so far

# Define dimensions of X 
N=700;
T=900;
lambda=N/T;

# Randomize X
x=randn(N,T);
# Estimate standard deviation
s=std(x(:));
# Now we estimate y and it's eigenvalues
y=x*x'/T;
l=eig(y);
# Calculating lambda_plus and lambda_minus
lambda_plus=(s^2)*(1+sqrt(lambda))^2;
lambda_minus=(s^2)*(1-sqrt(lambda))^2;
x_lim = linspace(0.1, 3, 10^3)
# Taking into account indicator function
  if (x_lim <= lambda_plus && x_lim>= lambda_minus);
    smth=(1./(2*pi*lambda*x_lim*s^(2))).*sqrt((lambda_plus-x_lim).*(x_lim-lambda_minus));
  else 
    smth = 0
  endif
  # Normalize histogram
  [f, p] = hist(l, nbins = 20)
 hold on;
   bar(p, f / sum(f));
 plot(x_lim, smth)

The problem I have with this function is that the output looks like following:

which doesn't exactly look as should. According to wikipedia these two plots should converge much more significantly. I found somewhere on internet very similar plots and those look like this:

My question is: what I'm doing wrong ? Did I choose a wrong scale ? As I can see on the second plot function takes much less values than mine but I have no idea what is wrong with its calculation. Could you please give me a hand telling what I did incorrectly ?