Problem in theoratical and simulated bit error rate

Dear fellows i am facing one problem.Problem is that there is difference between simulated and theoretical curve while calculating BER.I mean when i plot both curves both of them are separated but actually both of them should overlap each other.i am pasting my code here.

clc,close all,clear all
codn=1000;
EbN0=30:50;
q=1.6e-19; %
% signal-to-noise ratio in dB.

Ib=202e-6; % Background Noise Current+interfernce
N0=2*q*Ib; % Noise Spectral Density, 2*q*Ib
fs=360;
bode=10;
Resp = 1; %Receiver responsivity 
bitRate = 10^1 %BitRate 
Tbit = 1/bitRate; %one bitTime 
code=round(rand(1,codn))
code_len=round(1/bode/(1/fs)) % no of samples/symbol
SNR=10.^(EbN0./code_len);
sgma=zeros(1,length(SNR));
ber=zeros(1,length(SNR));
for i=1:length(SNR)
    Pavg(i) = sqrt((N0*bitRate*SNR(i))/(2*Resp^2)); %Luminous power     
    Ipeak(i) = 2*Resp*Pavg(i); %Photodiode Current    
    Epeak(i) = Ipeak(i)^2 * Tbit; %Peak current energy     
    sigma(i)=sqrt(N0*Epeak(i)/2); %standard deviation after receiver    
    threshold=0.5*Epeak(i); %threshold level 
 for ii=1:codn  
     x((ii-1)*code_len+1:code_len*ii)=code(ii);
 end
 cl=SNR(i)
 size(x)
  
%  y=x+sgma(iii)*randn(1,length(x)); % if want to add noise
 x2 = x-(1/2);    % get rid of most of the dc peak
pt=ones(1,code_len);
rt=pt;
% set up time and frequency arrays
length(x)
 u = length(x2);
 N = 2^ceil(log2(abs(u)));
 
delt = 1/fs;
delf1=fs/u;
figure(1)
tvec2=(1:length(x2))*delt;
plot(tvec2,((x2(1,:)+0.5)))
title('orignal baseband')
 xlabel('time');
 ylabel('amplitude')
ylim([-1 1.5]);
 y = fftshift(fft(x2)/N);  
z=abs(y);
figure(2)
fvec2=(-length(x2)/2:length(x2)/2-1)*delf1;
plot(fvec2,z)
title('FFT')
xlabel('frequency')
 ylabel('amplitude')
 figure(3)
 z=y;
 z(abs(fvec2)>50 & abs(fvec2)<=150)=0;
 plot(fvec2,abs(z))
 xlabel('frequency removed from 50 to 150 HZ');
ylabel('amplitude')
figure(4)
zf=fftshift(z)*N;
zifft=ifft(zf)+0.5;
MF_out=conv(zifft,rt)*0.5; % sampling time = 0.1
MF_out_downsamp=MF_out(code_len:code_len:end);
MF_out_downsamp=MF_out_downsamp(1:10);
Rx_th=zeros(1,codn);
%Rx_th(find(MF_out_downsamp>1/2))=1;
%Rx_th(find(MF_out_downsamp>mean(abs(MF_out_downsamp))))
   for k=1:length(codn)
    receivedSignal(k) = code(k)*Epeak(i)+ normrnd(0,sigma(i));
    end

                    if MF_out_downsamp(k)>mean(abs(MF_out_downsamp))
    Rx_th(k)=1;
                end
                    if code(k)>mean(abs(code))
                        code(k)=1;


end
nra=code==Rx_th
nerr=length(nra)-sum(nra)
ber(i)=1/10^(nerr/length(nra)*cl)
end
plot(tvec2,(abs(zifft)))
ylim([-1 1.5])
title('recovered signal')
xlabel('time');
 ylabel('amplitude')
 theorBER = qfunc(sqrt(SNR)); %theorical formula of OOK BER
semilogy(EbN0,theorBER,'g',EbN0,ber,'b') %theoretical BER graph 
legend('theory','theory');
grid on  
ylabel('BER') 
xlabel('SNR (dB)')
title('Bit Error Rate ')  %graph definition 


legend('simulation','theory');