Simulating a 2D- Gaussian field. Results are not what I expected

Hi guys. So basically i want to do somthing simulate a 2D gaussian field from a RBF kernal where i can specify the hyper paramters and and area for which the field is generated over. Basically I want to do exactly this

Went into this thinking it should be simple but either I was wrong or im missing somthing. What im getting is slighty off what i want

the first image is for a low length scale and the second is for a length scale of 100. I dont understand why its stretching on the Y-axis like that

N=200;
x=linspace(0,10,N);
y=linspace(0,10,N);

[X,Y]=meshgrid(x,y);

%covarience matrix from RBF kernal
covm=Dkernal(1,X,Y,100);
%fix pos semi def error
covm = covm+.0001 * eye(N);

%generate gaussian process
R=chol(covm);
z = randn(N,N);
mu = zeros(N,1)';
x = mu + R'*z;

imagesc(x);

function [cov] = Dkernal(sigma,x,y,l)
    cov=sigma^2*exp(-(squareform(pdist(x.'))).^2/(2*l^2))+exp(-(squareform(pdist(y.'))).^2/(2*l^2));
end

Thats the code im using. Im pretty sure the error is somwhere in the kernal function but i cant be sure. Any advice on this would be appreciated