clear;

[X,Y] = meshgrid(-100:1:100,-100:1:100);
x1 = 20; y1 = 20; x2 = -50; y2 = 10; sig1 = 15; sig2 = 15; C1 = [sig1^2 0; 0 sig1^2]; C2 = [sig2^2 0; 0 sig2^2];
Z = 0.5*exp(-((X-x1).^2+(Y-y1).^2)/(2*sig1*sig1)) + 0.5*exp(-((X-x2).^2+(Y-y2).^2)/(2*sig2*sig2));
figure, surf(X,Y,Z); shading flat;

N = 2000;
a = zeros(2000,2);

for i=1:N
    v = rand(1); 
    if v <= 0.5
        a(i,:) = mvnrnd([x1 y1],C1);    % samples from a 2-dimensional Gaussian with mean vector [x1 y1] and covariance matrix C1 defined above
    else
        a(i,:) = mvnrnd([x2 y2],C2);    %  samples from a 2-dimensional Gaussian with mean vector [x2 y2] and covariance matrix C2 defined above       
    end
end

figure, scatter(a(:,1),a(:,2));