% t = linspace(0, 20*pi, 1000);
% c = cos(t);
% s = sin(t);
% plot3(c, s, t)

% xs=-3:0.01:3;
% mu = 1;
% sigma = 1;
% plot(xs,normpdf(xs,mu,sigma));


% t = [0 : pi/100 : 20*pi];
% c = cos(t);
% s = sin(t);
% plot3(t.*c, t.*s, t)

% x = [-3:0.1:3]';
% y = [-2:0.1:2]';
% [X, Y] = meshgrid(x, y);
% F = (X + Y).*exp( -X.*X - 2*Y.*Y );
% mesh(X, Y, F)
% surf(X, Y, F)

mu    = [1 1]';   
sigma = [3 4; 4 3];
S = inv(sigma);
% 
% x = [0:0.2:6]';
% y = [0:0.2:6]';
% [X, Y] = meshgrid(x, y);
% 
% X = X - mu(1);
% Y = Y - mu(2);         
% 
%(x y) (s1 s2) ( x
%       s3 s4)   y)
% = ( x*s1+y*s3, x*s2+y*s4)*(x y)'
% =   x*x*s1+x*y*x3 + x*y*s2+y*y*s4
% 
% Z=(1/(2*pi*abs(det(S))^(1/2)))*exp((-1/2)*(X^2*S(1,1)+X*Y*S(2,1)+X*Y*S(1,2)+Y^2*S(2,2)));
% mesh(x, y, Z)
% surf(x, y, Z)

x=0:9; y=0:9;
[X, Y]=meshgrid(x,y)

i=1:10; j=1:10;
for i=1:10
    for j=1:10;
z=[x(i) y(j)];
f=z*inv(sigma)*z';
(1/(2*pi*abs(det(S))^(1/2)))*exp((-1/2)*f)
end
end