sigmax = 3;
sigmay = 2;
rhoxy = 0.5;

tmp = 8;
point_n = 31;
x = linspace(-tmp, tmp, point_n);
y = linspace(-tmp, tmp, point_n);
[xx, yy] = meshgrid(x, y);

tmp = (xx/sigmax).^2 - 2*rhoxy*(xx.*yy)/(sigmax*sigmay)+(yy/sigmay).^2;
expo = tmp*(-1/(2*(1-rhoxy^2)));
zz = exp(expo)/(2*pi*sigmax*sigmay*sqrt(1-rhoxy^2));
subplot(2,2,1);
mesh(xx, yy, zz);
axis([-inf inf -inf inf -inf inf]);
set(gca, 'box', 'on');
%frot3d on
rotate3d on

subplot(2,2,2);
contour(xx, yy, zz, 15);
axis image;