Index in Position 1 exceeds array bounds (must not exceed 100).

alpha = 0.5;
 w = gamma_rdn(M,alpha);
 x1 = (0.0001:0.001:1); % For plot
 
 figure(5)
 subplot(2,1,1);hist(w);title('Histogram of Gamma RDN');
 subplot(2,1,2);plot(x1,pdf('gam',x1,alpha,1));title('Theoretical Gamma Density with \alpha = 0.5');
 axis([0 1 0 100]);
 
 % The gamma_rdn function is implemented as follows:
 function[w] = gamma_rdn(M,alpha)
  % Generate random numbers from the gamma distribution with parameter
  % alpha <= 1, beta = 1
  pe = exp(1);
  w = zeros(M,1);
  u = rand(100,1);
  b = (alpha + pe)/pe;
  i = 0;
  j = 0;
  while j < M
  i = i+1;
  y = b*u(i,1);
  if y <= 1
  z = y^(1/alpha);
  i = i+1;
  if u(i,1) <= exp(-z)
  j = j+1;
  w(j,1) = z;
  else
  i = i+1;
  end
  else
  z = -log((b-y)/alpha);
  i = i+1;
  if u(i,1) <= z^(alpha - 1)
  j = j+1;
  w(j,1) = z;
  else
  i = i+1;
  end
  end
  end
  if i > 95
  u = rand(100,1);
  i = 0;
  end
 end