% The number of samples we are going to use.
n = [100 ; 1000; 10000];
% Loop over each
for idx = 1:3
  % Generate the random deviates
  y = poissrnd(5*ones(n(idx),1));

  % The max of y is how far out we need to go in terms of 
  % estimatng the PDF
  my = max(y);x = (0:my)';
  % For loop to build up the reltive frequency information
  ny = zeros(my+1,1);
  for nidx = 0:my
    ny(nidx+1) = sum(y==nidx)/n(idx);
  end
  % Generate the "theoretical" values for the PDF
  p = poisspdf(x,5);
  % Plot and print the results
  plot(x,p,'-',x,ny,'--');
  xlabel('n')
  title(['n = ',num2str(n(idx))]);
  legend('Theory','Sample');
  pstr = ['print -djpeg p2p10p5_',num2str(idx),'.jpg'];
  eval(pstr);
end