Basics of AWGN and SNR signal simulation in MATLAB

I'm trying to get my head around signal simulation in MATLAB/octave and generation of AWGN using

1:

    n = sqrt(1/EbN0)*randn(1,N);

vs 2:

    n = sqrt(0.5 * (1/EbN0))*randn(1,N);

and

3:

    n = sqrt(1/EbN0(i))*(randn(1,N) +1i*randn(1,N));

where `EbN0` is my Bit Energy to Noise PSD ratio.

In each of these cases, I want to keep my Signal Power as 1 and change Noise Power.

I understand that `rand()` gives values with variance as `1`.

So, when we want both real and imaginary part, we use:

    n = sqrt(1/2)*(randn(1,N) +1i*randn(1,N));

To bring the variance down to 1 from 2. But how does this scale to `EbNo` values?

I mean, going by this logic, shouldn't I be dividing my $\sqrt{\frac{N_0}{2}}$ to get total Power to $N_0$?

If that is the case, then, as $N_0 = \frac{1}{E_bN_0}$, we should be multiplying the normalized noise by $\frac{1}{2B_bN_0}$?

So, how does it work?

When I use the relation in (1), my simulations tend to disagree with theoretical BER-SNR plots (even when the variance of real-valued noise is 1). If I use (2), the plots match.

Consider the code:

    N = 10^5;
    x = randi([0,1],1,N);
    x = (x*2)-1;

    SNR_dB = -10:1:15;
    EsN0 = 10.^(SNR_dB/10);

    SER = zeros(1,length(SNR_dB));
    for i = 1:length(SNR_dB)
        n = sqrt(0.5 * (1/EsN0(i)))*randn(1,N);
        y = x + n;
        SER(i) = length(find(y.*x<0))/N;
    end

    SER_Theory_Q = qfunc(sqrt(2*EsN0));
    semilogy(SNR_dB, SER_Theory_Q, 'g--s', 'Linewidth', 2);
    hold on;

    SER_Theory_E = 0.5 * erfc(sqrt(EsN0));
    semilogy(SNR_dB, SER_Theory_E, 'm--o');

    semilogy(SNR_dB, SER, 'k->');

Here, I'm using the variant (2) and the results agree closely with the theoretical values.

If I use (1) instead, the plots don't agree.