Sigma = rotMat([1,1])*[3 0;0 1]*rotMat([1,1])'; c = 1; n = 101; x = linspace(-3,3,n); y = linspace(-3,3,n); [X,Y] = meshgrid(x,y); Z = zeros(n,n); fx = zeros(n,n); for i = 1:n; for j = 1:n Xnow = [X(i,j); Y(i,j)]; Z(i,j) = sqrt(Xnow'*inv(Sigma)*Xnow); fx(i,j) = exp(Xnow'*inv(Sigma)*Xnow*(-1/2))/sqrt(2*pi*det(Sigma)); end end figure(1);clf; surf(X,Y,fx); xlabel('x_1'); ylabel('x_2'); zlabel('f(x)') view([-5,42]) print -dpng mvn1 figure(2);clf; [C h]=contour(X,Y,-Z); set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2) xlabel('x_1'); ylabel('x_2'); print -dpng mvn2 hold on; plot(sqrt(3)*[0 1/sqrt(2)], sqrt(3)*[0 1/sqrt(2)],'r','LineWidth',2); plot(sqrt(1)*[0 -1/sqrt(2)], sqrt(1)*[0 1/sqrt(2)],'b','LineWidth',2); axis equal; print -dpng mvn3