r = rand(1000,2); % 2 cols of uniform rand
% Do the Box-Muller transformation
n = sqrt(-2*log(r,1)))*[1,1].*[cos(2*pi*r,2)), sin(2*pi*r,2))];
hist % plot two histograms
code 2:
function x=boxmuller(n,sd);
x=zeros(n,1);
len=ceil(n/2);
rand('seed',sd);
for i=1:len
u1=rand(1,1);
u2=rand(1,1);
x(i)=sqrt(-2*log(u1))*cos(2*pi*u2);
if (i+len<=n)
x(i+len)=sqrt(-2*log(u1))*sin(2*pi*u2);
end;
end;