Need help for AFsimulation in AWGN channel

Hi, everyone, I want to compute the BER/SER of the AF relay protocol in AWGN channel with QPSK. But i don't get the rigth result. the codes are the following:

mode=4; % QPSK;
SNRdB =0:2:40; % SNR(dB)
NumBits=1e6; % number of total bits
BPS=log2(mode); % bit per symbol
data=randint(1,NumBits); % datas
Txsig=pskmod(data,mode); % % QPSK modulation
%===========================================
for i =1:length(SNRdB)

SNR=10.^(SNRdB(i)/10);
% CN-(0,1)
N_sd=(randn(1,NumBits)+j*randn(1,NumBits))/sqrt(2)*10.^(-SNRdB(i)/20);
N_sr=(randn(1,NumBits)+j*randn(1,NumBits))/sqrt(2)*10.^(-SNRdB(i)/20);
N_rd=(randn(1,NumBits)+j*randn(1,NumBits))/sqrt(2)*10.^(-SNRdB(i)/20);
% %==================================================

% % AWGN channel is it right???
chan_sd=awgn(Txsig,SNRdB(i)); % s-d channel
chan_sr=awgn(Txsig,SNRdB(i)); % s-r channel
chan_rd=awgn(Txsig,SNRdB(i)); % r-d channel

Ysr=chan_sr.*Txsig+ N_sr; %
belta=sqrt(SNR)./(sqrt(SNR.*abs(chan_sr).^2)+1);
Y=Ysr.*belta; %
Ysd=chan_sd.*Txsig+ N_sd;
Yrd=chan_rd.*Y+ N_rd; %

% MRC
a1=conj(chan_sd);
a2=(conj(chan_sr).*belta.*conj(chan_rd))./(belta.^2.*abs(chan_sr).^2+1);
MRC=Ysd.*a1+Yrd.*a2; %

% QPSK demodulation
Rx_sd=pskdemod(Ysd,mode);
Rx_Mrc=pskdemod(MRC,mode);

% BER
[bitsd,bersd(i)]=biterr(data,Rx_sd); % compute the BER of the S-D link
[bitrd,bermrc(i)]=biterr(data,Rx_Mrc); % compute the BER of the S-R-D link
end
% ======================

% Pe1=berfading(SNRdB,'psk',mode,1); % compute the theoretical BER of the S-D link
% Pe2=berfading(SNRdB,'psk',mode,2); % compute the theoretical BER of the S-R-D link

%%=========================================
I want to know is the code right? After the simulation, the simulation curve can not coincide with the theoretical curve .But I can not find where the problem is. Please give me some help. Thank u!!