Generating non-homogeneous spatial Gaussian data I want to generate a spatial data following multivariate Gaussian distribution. 
However, I don't want it to be homogeneous, meaning I don't want the correlation/covariance to be homogeneous. I want it to be heterogeneous. 
Any suggestions on how to generate such data?
Here is an example of what I mean by heterogeneous covariance matrix
   25.8621   25.1207   24.6305   24.3867
   25.1207   25.3719   24.8768   24.6305
   24.6305   24.8768   25.3719   25.1207
   24.3867   24.6305   25.1207   25.8621

There are four variables. As you can see the marginal variance is different which means C(0) is dependent on location. So does C(h) for any h. This is one method of generating heterogeneous data.
But I want some spatial data consisting of like land, lakes which have very less covariance between each other.
 A: You might find this useful:
http://www.mathworks.com/help/stats/gmdistribution.fit.html
MU1 = [1 2];
SIGMA1 = [2 0; 0 .5];
MU2 = [-3 -5];
SIGMA2 = [1 0; 0 1];
X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];

scatter(X(:,1),X(:,2),10,'.')
hold on
options = statset('Display','final');
obj = gmdistribution.fit(X,2,'Options',options);
h = ezcontour(@(x,y)pdf(obj,[x y]),[-8 6],[-8 6]);

Now these covariances are only diagonal.  You can fix that by changing numbers.
The "gmdistribution.fit" is where you get the contour values.
So here is a plot that this would create:

Now you can see that it creates two multivariate gaussian distributions.  You only need to create one.
The function that generates the data is:
mvnrnd(MU1,SIGMA1,1000)

the mvnrnd, or mv-n-rnd is the multivariate normal random number generator.  Its inputs are the multivariate mean, covariance matrix, and desired sample count.  The output is an array of numbers of dimension informed by mu1.
The example covariance that you provided was 5 dimensional.  Here is something "hackable" that you should be able to convert for your own non-nefarious purposes.
%this creates a 5-dimensional mean
MU1 = rand(1,5)+[1,2,3,4,5];  

%this creates a 5x5 covariance, diagonally dominant
SIGMA1 = 23*rand(5,5)+eye(5);  

%this is for consistent notation.  previously it stacked multiple multivariate gaussians.
X = [mvnrnd(MU1,SIGMA1,1000)];  

%this plots the dots
scatter(X(:,1),X(:,2),10,'.');  

%this lets multiple plots overlay
hold on;  

%this is a parm of the gm fit function.  It says "display final result".
options = statset('Display','final'); 

%this fits a gaussian to the data
obj = gmdistribution.fit(X,1,'Options',options); 

%this plots the contours of the fitted gaussian over the domain of the data
h = ezcontour(@(x,y)pdf(obj,[x y]),[min(X(:,1)) max(X(:,1))],...
                                   [min(X(:,2)) max(X(:,2))]); 

Best of luck.
