Very nice and decent code for displaying scatterplot in the 2D map with color for empirical distribution. The counting part was efficient only if you have a larger dataset but set a smaller # of bins. If you increase bin to larger than certain value it starts to get slowly even you have smaller dataset. I can't think of anything better but use an alternative approach for the condition with smaller dataset, and switch to the original method if the # of bin is actually smaller.

Following is my approach, similar to BLu's but with better decision switching between two methods. BTW, the centering bug was fixed in this version based on Thomas's post so length(C) == n

% do counts
if numel(x) < n^2
% New method
binIntX = diff(limitX)/(n-1);
binIntY = diff(limitX)/(n-1);
for idx = 1:numel(x)
idxX = min([round(x(idx)/binIntX)+1 numX]);
idxY = min([round(y(idx)/binIntY)+1 numY]);
C(idxY,idxX) = C(idxY,idxX) + 1;
end
else
% Old method
for i = 1:numY-1
for j = 1:numX-1
C(i,j) = length(find(x >= xEdges(j) & x < xEdges(j+1) &...
y >= yEdges(i) & y < yEdges(i+1)));
end
end
end