There are multiple ways. histrogram2() would be the best alternative but it wasn't introduced until r2015b. I recreated the image you attached by using heatmap() but it could also be done using imagesc() or other methods. Look out for some small difference such as how the bins are labeled.
After loading the data I remove NaNs and then bin the data using histcounts(). Unfortunately heatmap() cannot label the edges of the bins; it only labels their centers. So I assign the data to the bin centers. The heatmap() function does the rest.
load carbig
% Remove NaNs
nanIdx = isnan(MPG) | isnan(Weight);
MPG(nanIdx) = [];
Weight(nanIdx) = [];
% bin the data
nBins = 10; %number of bins (there will be nBins + 1 edges)
[~, mpgEdge, mpgBin] = histcounts(MPG, nBins);
[~, wgtEdge, wgtBin] = histcounts(Weight, nBins);
% Calculate center of each bin
mpgBinCnt = mpgEdge(2:end) - (diff(mpgEdge)/2);
wgtBinCnt = wgtEdge(2:end) - (diff(wgtEdge)/2);
% Assign data to bins
MPGb = mpgBinCnt(mpgBin)';
Wightb = wgtBinCnt(wgtBin)';
% Put into table and plot heatmap
tbl = table(MPGb, Wightb);
hh = heatmap(tbl, 'MPGb', 'Wightb', 'CellLabelColor', 'none');
% Revers y axis
hh.YDisplayData = flipud(hh.YDisplayData);
% Label the axes
hh.XLabel = 'MPG';
hh.YLabel = 'Weight';
hh.Title = 'hist3 simulation';
% Change colormap
colormap parula
