the analytic density, the asymptotic approximation and the numerical computer-generated approximation for a given exponential distribution

expdist.M

N= 5000;
n = 2;
lambda = 1;
mu = 3;
p = ones([n,1]);
x = zeros(n,1);
e = zeros(n,1);
y = zeros(n,1);
for i =1:N;
   x(:,i) = exprnd(lambda,n,1);
   e(:,i) = x(:,i) - lambda*p;
   y(:,i) = mu*p + e(:,i);
   munhat = mean(y); 
   wn = sqrt(n)*(munhat - mu);
end

anlytcwn = inv(gamma(n)*lambda.^n)*sqrt(n)*((sqrt(n)*(wn + sqrt(n)*lambda)).^(n-1)).*(exp(-(sqrt(n)/lambda)*(wn + sqrt(n)*lambda)));
asympwn =  normpdf(wn,0, lambda.^2);
draws = wn;
%draws = anlytcwn;
%draws = asympwn;
gridsze = 300;
lowbd = min(draws); 
upbd = max(draws);
stderr = std(draws);
interquart = prctile(draws,75) - prctile(draws,25);
A = min(stderr,(interquart/1.34));
h = .9*A*(length(draws)^(-1/5));

dom = linspace(lowbd,upbd,gridsze)';
ran = zeros(gridsze,1); 

for i = 1:gridsze;
	a = -.5*sign( abs( (draws - dom(i,1))/h ) -sqrt(5) ) + .5;
	draws2 = nonzeros( draws.*a );
	tempp = zeros(length(draws2),1);	
		for j = 1:length(draws2);	
		tempp(j,1) = (.75/sqrt(5))*(1 -.2*( (draws2(j) - dom(i,1))/h )^2);
		end;	
ran(i,1) = (1/( (length(draws)*h) ) )*sum(tempp);
end;


plot(dom, ran)
hold on;
plot(wn, asympwn)
hold on;
plot(wn, anlytcwn)

the code is in Octave