Let's take this small example in which we generate just 3 pairs of correlated random numbers:
n=3; % Generate 3 pairs, instead of 1000
sigma=0.5;
Sigmalnd=sigma.^2.*[1 0; 0 1];
rng default % I corrected your syntax here
ZInd=mvnrnd([0 0],Sigmalnd,n)
The covariance matrix specifies that the first column and the second column will be uncorrelated. (Specifically, the parent distribution variables will be uncorrelated. The sample distribution will have some correlation, due to sampling error.)
Now let's look at the mean of that first row, as you suggested:
mean(ZInd(1,:))
You are saying that you want the mean of the first row to be zero. This implies that you want the sum of the first row to be zero. This further implies that the second column value is always the exact negative of the first column value. So, these values are not uncorrelated, but rather are perfectly (and negatively) correlated.
So, you just need to change Sigmalnd to sigma.^2.*[1 -1; -1 1]:
Sigmalnd=sigma.^2.*[1 -1; -1 1];
rng default
ZInd=mvnrnd([0 0],Sigmalnd,n)
Or, more simply, just generate the first column, and then negate it to get the second column.
rng default
ZInd=sigma*randn(n,1);
ZInd(:,2) = -ZInd(:,1)
