function X = ar_rnd(m,n)

% Implements a random number generator
% using the acceptance-rejection method.

f = inline('x.*exp(-(x.^2)/2)');  % pdf from which we wish to generate random numbers.
g = inline('exp(-x)');  % Comparison pdf: exponential distribution with mu = 1.

h = inline('x.*exp(x-(x.^2)/2)'); % h(x) = f(x)/g(x).
a = (1 + sqrt(5))/2; % h has a maximum at x = a.
c = h(a);  % Minimum value of c for which  c >= f(x)/g(x)  for all x.
                
X = zeros(m, n);

for i = 1:m*n
    accept = 0;
    while accept == 0
       r = exprnd(1);
       u = rand;
       if c*u <= f(r)/g(r)
           X(i) = r;
           accept = 1;
       end
    end
end