m = 31; n=m*m;   % number of grid-points, m > 40 are computationally demanding!
L = 4;                                  % limits of gridpoints are "-L" and "L"
noise = 0.001;                                              % variance of noise
[x1,x2] = meshgrid(-L:2*L/(m-1):L, -L:2*L/(m-1):L);
x1=reshape(x1,n,1); x2=reshape(x2,n,1);

Q=exp(-0.5*abs(x1*ones(1,n)-ones(n,1)*x1').^2 ...     % this is the covariance
      -0.5*abs(x2*ones(1,n)-ones(n,1)*x2').^2);         

y=randn(1,n)*chol(Q+noise*eye(n));

surf(reshape(x1,m,m),reshape(x2,m,m),reshape(y,m,m))

