How can I map a Gaussian vector according to an equation that has a spatial autocorrelation function?

I have the following function:

I want to run some Monte Carlo simulations and at each simulation I have to generate new values of E at certain locations x and z. So basically I will have to generate a grid of points with coordinates and where I will evaluate E.

is a 0 mean Gaussian random field with a variance equal to 1. It has a spatial autocorrelation function where is the horizontal distance between two points and is the vertical distance between two points. The function R has the folllwing equation:

I know how to generate the vector Y if it has no correlation but I'm not sure how to do it with the autocorrelation function, I did the following but I'm not sure if this is correct:

Generated a grid of points (pts), code below:

x=-50:1:50;
z=0:1:20;
pts=zeros(length(x)*length(z),2);

for i=1:length(x)
    for j=1:length(z)
        pts((length(z)*(i-1)+1)+j-1,:)=[x(i),z(j)];
    end
end

Generated a correlation matrix from the correlation function, code below:

cormat=zeros(length(pts),length(pts));

for i=1:length(pts)
        for j=i:length(pts)
            dx=pts(i,1)-pts(j,1);
            dy=pts(i,2)-pts(j,2);
            cormat(i,j)=exp(-abs(dx)/10.-abs(dy)/3.);
        end
end

Generated a covariance matrix from the correlation matrix:

cormat=cormat'+cormat-eye(length(pts));
stdmat=eye(length(pts));
covmat=stdmat*cormat*stdmat;

Evaluated the function Y using MATLAB function mvnrnd:

y=mvnrnd(ones(1,length(pts)),covmat);