Accept/Reject Method for Monte Carlo Simulation

I have a probability distribution function given by f(m) = summation function from i = 0 to infinity of a_i * P_i(m), where P_i represents Legendre polynomials, and m is equal to cos(theta) and ranges from -1 to 1. I am given the first four coefficients and polynomials, shown below.

a0 = 1, a1 = 1/2, a2 = 1/3, a3 = 1/4

P0(m) = 1, P1(m) = m, P2(m) = 1/2(3m^2 - 1), P3(m) = 1/2(5m^3 - 3m)

I am trying to generate 10^4 random numbers to represent the probability density function, using the rejection sampling method. The code I am using is given below.

clear all;
close all;

Numtrials = 1e3;                  % number of trials 
rn_attempts = zeros(Numtrials,1); % array of attempted random numbers 
NRN = 1e4;                        % number of random numbers 
RN_accept = zeros(NRN,1);         % array for accepted random numbers 
nbin = 100;                       % number of bins 
a = -1;                           % lower bound of target random number 
b = 1;                            % upper bound of target random number 

% positions of each bin, use the center of the bin 
xbin = zeros(nbin,1); 
bin_length = (b-a)/nbin;
analytical_pdf = zeros(nbin,1); % analytical pdf at the center of each bin

for i=1:nbin 
    xbin(i) = (i-0.5)*bin_length+a;
    analytical_pdf(i) = probdist(xbin(i));
end 

pdf_max = 25/12; %sum of all coefficients a_i 

n = 0; 
nacc = 0; % # of accepted random numbers 
while(n < Numtrials & nacc < NRN) 
    n=n+1; % track the trial times
    % generate a random number in (0 1) and linearly map it into [1 b]
    rand1 = rand(1);
    
    % map and store the random number 
    rand1 = a+rand1*(b-a);
    rn_attempts(n) = rand1;
    
    probab = probdist(rand1)/pdf_max; 
    
    % generate second random number to compare to scaled pdf 
    rand2 = rand(1);
    if(rand2 < probab) 
        nacc = nacc+1;
        
        % store the random number if accepted 
        RN_accept(nacc) = rand1;
    end 
end 
    
% counts of accepted and tried random numbers in each bin 
counts_try = hist(rn_attempts,nbin);
counts_accept = hist(RN_accept,nbin); 
rate_accept = counts_accept./counts_try; 


% plot result 
plot(xbin,rate_accept,'ro',xbin,analytical_pdf/pdf_max,'g-');
axis([-1 1 -1 1]);
title('Rejection sampling');
xlabel('x');
ylabel('Scaled probabilty distribution function');
legend('Rejection sampling','Analytical','Location','Northwest');
set(gca,'fontsize', 12);

function f=probdist(m) 
    a0 = 1;
    a1 = 1/2;
    a2 = 1/3;
    a3 = 1/4;
   f=a0+a1*m+a2*0.5*(3*m^2-1)+a3*0.5*(5*m^3-3*m);
end 

However, the graph I get shows the accepted rate to be less than and greater than the analytical solution. What could I do to correct this?