Simulation of an hyper geometric random variable

Function code:

 
function [y]=HyperGeom(m,n,Mean,var)
//------------------------------- 
z = var / (Mean * Mean);
pP = 0.5 * (1.0 - sqrt((z - 1.0) / ( z + 1.0 )));
y=rand(m,n,'uniform')
zz=find( y > pP) ;
y=pP*ones(y);
y(zz) = (1-pP)*ones(zz);
y1=rand(m,n,'uniform')
y=-Mean * log(y1) ./ (2.0 * y) ;

Simulation code:

n=10000; 
y=HyperGeom(1,n,1,10)
histplot([0:0.25:10],y,[1,1],'061',' ',[0,0,10,0.4]);