Steve Simon

This would be really useful for scatter3 too, any suggestions how I could do that?

It works well for small datasets and small n (such as 50).
I added a none option, as suggested above and changed the centering problem as suggested above.
I tried a large dataset and n=200, and have to terminate the calculation because it is taking too long. The code used in calculating the counts is not the most efficient. cloudPlot has a much more efficient algorithm. I made modification based on cloudPlot.

error(nargchk(2,6,nargin),'struct');

x = x(:);
y = y(:);

%new code
pointSelect = isinf(x) | isnan(x) | isinf(y) | isnan(y);
x = x(~pointSelect);
y = y(~pointSelect);
%new code ends

if length(x) ~= length(y)
    error('SCATTERCLOUDDataVectorSizesDoNotMatch','The number of elements in x and y do not match')
end

if nargin < 6
    %cmap = flipud(gray(256));
    cmap='default';
end

if nargin < 5
    %clm = 'k+';
    clm='none';
end

if nargin < 4
    l = 1;
end

if nargin < 3
    n = 25;
end

% min/max of x and y
minX = min(x);
maxX = max(x);
minY = min(y);
maxY = max(y);

% edge locations
xEdges = linspace(minX,maxX,n);
yEdges = linspace(minY,maxY,n);

% shift edges
xDiff = xEdges(2) - xEdges(1);
yDiff = yEdges(2) - yEdges(1);
xEdges = [-Inf, xEdges(2:end) - xDiff/2, Inf];
yEdges = [-Inf, yEdges(2:end) - yDiff/2, Inf];

% number of edges
numX = numel(xEdges);
numY = numel(yEdges);

%new code
xBinIndex = floor((x - minX+xDiff/2)/xDiff)+1;
yBinIndex = floor((y - minY+yDiff/2)/yDiff)+1;
%new code ends
    
% hold counts
C = zeros(numY,numX);

% do counts
%new code
for i=1:numel(x)
    C(yBinIndex(i),xBinIndex(i))=C(yBinIndex(i),xBinIndex(i))+1;
end
%new code ends

%for i = 1:numY-1
% for j = 1:numX-1
% C(i,j) = sum(x >= xEdges(j) & x < xEdges(j+1) &...
% y >= yEdges(i) & y < yEdges(i+1));
% end
% disp(num2str(i));
%end

% get rid of Infs from the edges
xEdges = [xEdges(2) - xDiff,xEdges(2:end-1), xEdges(end-1) + xDiff] + xDiff/2;
yEdges = [yEdges(2) - yDiff,yEdges(2:end-1), yEdges(end-1) + yDiff] + yDiff/2;

%xEdges = [xEdges(2) - xDiff,xEdges(2:end-1), xEdges(end-1) + xDiff];
%yEdges = [yEdges(2) - yDiff,yEdges(2:end-1), yEdges(end-1) + yDiff];

% smooth the density data, in both directions.
C = localSmooth(localSmooth(C,l)',l)';

% create the graphics
ax = newplot;
s = surf(xEdges,yEdges,zeros(numY,numX),C,...
         'EdgeColor','none',...
         'FaceColor','interp');
view(ax,2);
colormap(ax,cmap);
grid(ax,'off');
holdstate = get(ax,'NextPlot');
set(ax,'NextPlot','add');

if ~strcmpi(clm,'none')
    p = plot(x,y,clm);
end

axis(ax,'tight');
set(ax,'NextPlot',holdstate)

% outputs
if nargout
    h = [s;p];
end

function B = localSmooth(A,L)
r = size(A,1);
I = eye(r);
D1 = diff(I);
D2 = diff(I,2);
B = (I + L ^ 2 * D2' * D2 + 2 * L * D1' * D1) \ A;

To fix the glaring centering bug, change these lines from:

xEdges = [xEdges(2) - xDiff,xEdges(2:end-1), xEdges(end-1) + xDiff];
yEdges = [yEdges(2) - yDiff,yEdges(2:end-1), yEdges(end-1) + yDiff];

TO:

xEdges = [xEdges(2) - xDiff,xEdges(2:end-1), xEdges(end-1) + xDiff] + xDiff/2;
yEdges = [yEdges(2) - yDiff,yEdges(2:end-1), yEdges(end-1) + yDiff] + yDiff/2;

In the future, if there's a "simple" fix, why not post it and be helpful?