n   = 10000;
rho = 0.5;

z1  = normrnd(0, 1, n, 1);
z2  = normrnd(0, 1, n, 1);

x1  = z1;
x2  = rho*z1 + (1 - rho^2)^(1/2)*z2;

y1  = power(x1, 2);
y2  = power(x2, 2);
corr(y1, y2)