BPSK bpsk in awgn not matching theoretical results

I am doing a monte carlo simulation for BPSK signals and I am getting a result of (roughly):

I know the result is supposed to be:

for a bpsk channel (source: https://www.mathworks.com/help/comm/ug/bit-error-rate-ber.html )

Can someone please look at my code provide some pointers as to what I'm doing wrong?

Fs = 10;               %Sample Frequency
Ts= 1/Fs;
t= 0:Ts:1;
degree=1;               %Degree of the space ie. BPSK= 1 
sim_length= 40000;      %Number of symbols sent in the simulation

error_count=0;
       %Define the two symbols to transmit
       p{1,1}=      sqrt(2)*sin(2*pi*t);
       p{1,2}=      sqrt(2)*sin(2*pi*t+pi);
 
x=0;
for Eb_No= 0:10       
       
       %Begin Monte Carlo for the desired SNR
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
            for i=1:sim_length 
                
          %ceil rounds up always
            symbol= ceil(2*degree*rand);            %Pick rand Symbol index
            symbol_to_tx= cell2mat(p(symbol));      %Create that symbol to transmit 

  %%%%%%%%% https://www.mathworks.com/help/comm/ug/awgn-channel.html
            Eb_No_SAMPLE  =     Eb_No+10*log10(1/Fs) ;



            noisy_sig= awgn(symbol_to_tx, Eb_No_SAMPLE, 'measured' );        %Signal is injected with noise
                                                           
 
 for j=1:2*degree
        %%Multiply against each symbol in the dicationary, and then
        %%integrate over one period 
        symbol_dictinory= cell2mat(p(j));
        test_funct= noisy_sig .* symbol_dictinory;
        ML_vector(j)=trapz(test_funct)* (1/Fs);   
            % (1/Fs) is required to normalize because of sampling
            % note that trapz is an approximation

 %Highest output in ML vector is the "most likely" symbol that was tx'ed.

 end

[~,ML_symbol] = max(ML_vector);

    %If the ML decoder was wrong, incrment the error count
    if ML_symbol ~= symbol
        error_count=error_count+1;
    end
P_error= error_count/sim_length; 


            end
       %End Monte Carlo for the desired SNR
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


x=x+1;
ber(x)=P_error;



%https://www.mathworks.com/help/comm/ug/bit-error-rate-ber.html
%%%This is just the theoretical error of the AWGN QPSK to compare
 Eb_No_db(x)=Eb_No;
 Eb_No_lin=10.^(Eb_No_db./10);
%    theory_ber=0.5.*erfc(sqrt(Eb_No_lin));
  theory_ber= qfunc(  sqrt(2*Eb_No_lin)   );



%%%Reset counters for next dB iteration
error_count=0;
P_error=0;

end



semilogy(Eb_No_db,ber,'-bo',Eb_No_db,theory_ber,'-mh')
grid on
title('BPSK with AWGN');
xlabel('E_b/N_o');
ylabel('P_b = BER');
legend('Simulated BPSK','bpsk theory')