16 QAM simulation over Rayleigh channel

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.

      % Create Rayleigh fading channel object.
      chan = rayleighchan(1/1000,100);

      g = chan.PathGains                          % Retrieve the PathGains property
      %chan.NormalizePathGains=0;

        dop_flat = doppler.flat
        chan.DopplerSpectrum = dop_flat
       chan.StorePathGains = 1;
      %chan.StoreHistory = 1;                              % Allow states to be stored
      chan.ResetBeforeFiltering=0;
      chan

      % Generate data and apply fading channel.
      M = 16;                                                                     % QAM  modulation order
      k = log2(M);
      hMod=comm.RectangularQAMModulator(M);               % Create a 16 QAM modulator
      %hMod.InputType='bit';
      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
      %tx1=dec2binvec(tx,5*1e4);

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

      %scatterplot(qamSig)
       axis([-20 20 -20 20]); % Set axis ranges.
      grid on;

       %update(scatterPlot,qamSig,'b');

      % 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) 
          %clear hDemod;
          hDemod=comm.RectangularQAMDemodulator;               % Create a 16 QAM Demodulator
          fadedSig = filter(chan,qamSig);                 % Apply the channel effects
          clear snr;
      clear rxSig1;

          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(n)] = 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