Multivariate Guassian Distribution

You're really trying to do two things here. The first is, you have some random data and you want to fit it to a multivariate normal distribution. Your approach to this part works, although it can be streamlined:

n = 1000; d=2;
X = randn(n,2);

Get mean and covariance:

mumat=mean(X);
cov_mat=cov(X);

The second part is plotting the resulting distribution. Here you need a regular grid for your variables, not the random values you generated above:

x = -3:.2:3; y = -3:.2:3;
[X,Y] = meshgrid(x,y);
X = X-mumat(1); Y = Y-mumat(2);

Combine X and Y in a way that each row represents one 2D variable.

Z = [X(:) Y(:)];

Now calculate the probabilities.

scale=((2*pi)^(d/2))*sqrt(abs(det(cov_mat))); 
p = zeros(length(Z),1);
for ii=1:length(Z)
  p(ii) = exp(-Z(ii,:)/cov_mat*Z(ii,:)'/2)/scale;
end

Reshape and plot.

p = reshape(p,length(x),length(y));
surf(x,y,p)
xlabel('x'), ylabel('y')

For example 

Lets Compute and plot the pdf of a bivariate normal distribution with parameters mu = [0 0] and Sigma = [0.25 0.3; 0.3 1].

Define the parameters mu and Sigma.

mu = [0 0];
Sigma = [0.25 0.3; 0.3 1];%initialisation

Now Creating a grid of evenly spaced points in two-dimensional space

x1 = -3:0.2:3;
x2 = -3:0.2:3;
[X1,X2] = meshgrid(x1,x2);
X = [X1(:) X2(:)];
Now we will evaluate the pdf of the normal distribution at the grid points.

y = mvnpdf(X,mu,Sigma);
y = reshape(y,length(x2),length(x1));
Now plotting 

surf(x1,x2,y)
caxis([min(y(:))-0.5*range(y(:)),max(y(:))])
axis([-3 3 -3 3 0 0.4])
xlabel('x1')
ylabel('x2')
zlabel('Probability Density')