Creating a Z-matrix of counts with the same range as my 2-D scatter plot

I have two vectors, X and Y, that I am using to create a traditional, vanilla scatterplot. However, the plot is very dense—about 3.3M points total. Both axes and the elements in the vectors have limits of [-256 256].

I would like to create a matrix, Z, that is 512x512 and corresponds to the scatter plot. I want the number of points falling at i,j on the scatter plot to correspond to the integer value Z(i,j). This is all for the purpose of using functions like surf() and mesh() with this Z matrix (right now I'm getting memory errors because the vectors are very large).

Below is the code that I currently have written. It takes at least 20 minutes to run, and I'm hoping for some suggestions to optimize it. Please let me know if you have any suggestions!

load('error_mat.mat')  % 3.3M x 2 array
err_x = round(error_mat(:,1));  % round these puppies to integer values
err_y = round(error_mat(:,2));
goalmat = zeros(512,512,'single');  % Z-matrix

for i = -255:256
    for j = -255:256
        potential_match_x = find(err_x == i);
        potential_match_y = find(err_y == j);
        ct = numel(intersect(potential_match_x,potential_match_y));  % intersect b/c each row of [err_x, err_y] corresponds to 1 point
        goalmat(i+256,j+256) = ct;  % store the count value
        clearvars potential_match_x potential_match_y 
    end
end

[Image Analyst]

You forgot to attach your data 'error_mat.mat' or any sort of visualization at all. So about I can suggest is to try histogram2().

[ohmstead]

Just posted an example of the scatter plot and a small subset of the error matrix. And I think you misunderstood my question. I'm not trying to plot a scatterplot. I'm trying to create a matrix representing the scatter plot.