Compressive sensing: numerical generation of RIP matrices

The restricted isometry property (RIP) states that:

$$\begin{equation} (1-\delta_K)||x||_2^2 \le ||A x||_2^2 \le (1+\delta_K)||x||_2^2 \end{equation}$$ for any $K$-sparse vector $x$ of length $N$. The corresponding restricted isometry constant is $\delta_K$, $0 < \delta_K < 1$.

I would like to perform a numerical experiment based this expression to produce a RIP-compliant matrix $A$ (Matlab code is supplied below). The matrix is constructed with Gaussian-distributed RV's (zero mean, variance $1/\sqrt{M}$). Having specified $N$ and $K$, the number of measurements $M$ is obtained with $M \geq K \log(N/K)$. This is a hit-or-miss procedure, I let the code run until the isometry constant falls between zero and one. Is this the proper way to generate $A$? Am I missing anything? Further analysis on $A$ (eigenvalues, etc.) doesn't seem to behave properly...

function [delta] = RIP_numerical(n,s)

% delta: restricted isometry constant; n: # columns; s: sparsity

% choose m (# measurements/rows) based on measurement criterion

m = ceil(s*log(n/s))


flag = 0; incr = 0; % initialize computing parameters

while flag ~= 1  % run until condition met

    incr = incr + 1
    % generate random, m x n normalized matrix N(0,1/m)
    A = (1/sqrt(m))*randn(m,n);
    sumA = sum(A);
    for k = 1:n
        A(:,k) = A(:,k)/sum(A(:,k));

    end

    % generate vector of length n, sparsity k
    x = zeros(n,1);
    list = randperm(n);
    for l = 1:s
        x(list(l))= randn;
    end

    y = A*x;    % obtain measurements


    norm_x = norm(x,2);  % l_2 norms
    norm_y = norm(y,2);

    delta = 1-(norm_y^2/norm_x^2);  % look at lower bound

    if (delta >0 & delta < 1) % condition met; save vars and kick out of routine.
        flag = 1
        save('RIP_config.mat', 'A','x','n','m','s');
    end

end