computing the covariance of a distribution

I am using a matrix to approximately represent a bivariate normal distribution. Here is the code:

x=[0:.1:10]; 
y=[0:.1:10];
[X,Y] = meshgrid(x,y);
prior = mvnpdf([X(:) Y(:)],[5,5],[1 0;0 1]); 
prior = reshape(prior,length(w1),length(w2));

As a result each element in the matrix represents a point on the distribution, each element's indices being its Cartesian coordinates.

How can I calculate the covariance matrix from this matrix. The function cov does not seem to be applicable, since it assumes that each column in the matrix is a variable, and that clearly is not the case with my matrix.

Please let me know if any clarification is needed and thank you for any help!