Simulation of an Erlang random variable

Function code:

 
function [y]=Erlang(m,n,pMean,pVariance) 
//------------------------------- 
  k = int( (pMean * pMean ) / pVariance + 0.5 );
  if (k <= 0) then k=1;end
  a = k / pMean;
 // we use blocks of size 100 to avoid overflows 
  res=[];
  ntir=100;
  ntirc=ntir;
  y=rand(ntir,k,'uniform');
  y= -log(prod(y,'r'))/a;
  res=[res;y];
  while ( ntirc < m*n )
     y=rand(ntir,k,'uniform');
     y= -log(prod(y,'r'))/a;
     res=[res;y];
     ntirc=ntirc+ntir;
  end 
  y=matrix(res(1:m*n),m,n);

Simulation code:

n=10000;
y=Erlang(1,n,10,1)
histplot(20,y,[1,1],'061');