16 QAM simulation over Rayleigh channel in Matlab

Hello everyone! I would like to demodulate 16 QAM modulated signal after Flat Rayleigh Channel and AWGN by usining modem.qamdemod. My BER graph is far from the thoretical. Any help will be appreciated.
Den.

Here is my code:

% Create Rayleigh fading channel object.
chan = rayleighchan;

g = chan.PathGains % Retrieve the PathGains property

dop_flat = doppler.flat
chan.DopplerSpectrum = dop_flat
chan.StorePathGains = 1;
chan.ResetBeforeFiltering=0;
chan

M = 16; % QAM modulation order
k = log2(M);
hMod=comm.RectangularQAMModulator(M); % Create a 16 QAM modulator
display(hMod)

%hDemod=comm.RectangularQAMDemodulator; % Create a 16 QAM Demodulator
%hDemod.outputtype='bit';
%display(hDemod)

% Initialize scatter plot.
scatterPlot = commscope.ScatterPlot;
scatterPlot.Constellation=16;

tx = randint(1e4,1,M); % Generate a random bit streams

qamSig = step(hMod, tx); % Modulate16 QAM the signal


% Compute error rate for different values of EbNo
EBNO = 1:50; % Range of EbNo values, in dB.
numEBNO = length(EBNO);


for n = 1:length(EBNO)

hDemod=comm.RectangularQAMDemodulator; % Create a 16 QAM Demodulator
fadedSig = filter(chan,qamSig); % Apply the channel effects

snr = n + 10*log10(k) - 10*log10(4);
rxSig1= awgn(fadedSig,n,'measured');

eqSig = rxSig1 .* conj(chan.PathGains);

rx = step(hDemod, eqSig); % Demodulate

% Compute error rate.
[nErrors, BER] = biterr(tx,rx);
BER

% Plot the new data
update(scatterPlot,eqSig,'r');

end;

% Compute theoretical performance results, for comparison.
BERtheory =berfading(EBNO,'qam',M,1)

% Plot BER results.
semilogy(EBNO,BERtheory,'b-',EBNO,BER,'r*-');
legend('Theoretical BER','Simulation BER');
xlabel('EbNo (dB)'); ylabel('BER');
title('16 QAM over Rayleigh Channel');
grid on;

end