%  2DGaussian

mu_x = input('Enter mean for X: ');
mu_y = input('Enter mean for Y: ');
sigma_x = input('Enter standard deviation for X: ');
sigma_y = input('Enter standard deviation for Y: ');

echo on

range = 4*max(sigma_x, sigma_y);
x = linspace(mu_x-range,mu_x+range,200)';
X = x*ones(1,200);
y = linspace(mu_y-range,mu_y+range,200)';
Y = ones(200,1)*y';

f_XY = 1/(2*pi*(sigma_x*sigma_y))*exp(-0.5*((X-mu_x).^2/sigma_x.^2 + (Y-mu_y).^2/sigma_y.^2));

mesh(X,Y,f_XY);
axis square
xlabel('X');
ylabel('Y');

pause

echo off
mu_x = input('Enter mean for X: ');
mu_y = input('Enter mean for Y: ');
sigma_x = input('Enter standard deviation for X: ');
sigma_y = input('Enter standard deviation for Y: ');
rho = input('Enter the correlation coefficient, rho: ');

range_x = 4*sigma_x;
range_y = 4*sigma_y;
x = linspace(mu_x-range_x,mu_x+range_x,200)';
X = x*ones(1,200);
y = linspace(mu_y-range_y,mu_y+range_y,200)';
Y = ones(200,1)*y';

echo on 

K = [sigma_x.^2 rho*sigma_x*sigma_y; rho*sigma_y*sigma_x sigma_y.^2];
f_XY = 1/(2*pi*sqrt(det(K)))*exp(-0.5*sum((([X(:)-mu_x, Y(:)-mu_y] * inv(K)) .* [X(:)-mu_x, Y(:)-mu_y])'));

f_XY = reshape(f_XY,200,200);

mesh(X,Y,f_XY);
axis square
xlabel('X');
ylabel('Y');

echo off