Optimize Reorganization of Matrix without Loops

TL;DR: I am trying to optimize the following short code in Matlab. Because it involves loops over large matrices, it is too slow.

for i = 1:sz,
    for j = 1:sz,
        if X(j) == Q(i) && Y(j) == R(i),
            S(i) = Z(j);
            break
        end
    end
end

Specifics: Basically, I started with three vectors of x, y and z data that I wanted to plot as a surface. I generated a mesh of the x and y data and then made a matrix for the corresponding z values using

[X, Y] = meshgrid(x, y);
Z = griddata(x, y, z, X, Y);

Because the data is collected in random order, when generating the surface plot the connections are all wrong and the plot looks all triangulated like the following example.

So, to make sure Matlab was connecting the right dots, I then reorganized the X and Y matrices using

[R, R_indx] = sort(Y);
[Q, Q_indx] = sort(X, 2);

From here I thought it would be a simple problem of reorganizing the matrix Z based on the indices of the sorting for matrix X and Y. But I run into trouble because no matter how I use the indices, I cannot produce the correct matrix. For example, I tried

S = Z(R_indx); % to match the rearrangement of Y
S = S(Q_indx); % to match the rearrangement of X

and I got this barcode...

Running the first block of code gives me the "desired" result pictured here. However, this takes far too long as it is a double loop over a very large matrix.

Question: How can I optimize this rearrangement of the matrix Z without for loops?

Please have a look at the following solutions, and test both with your matrices. Do they perform faster? The array indexing solution does, what you asked for, i.e. the re-arrangement of the matrices. The vector indexing might be even better, since it sorts your original vectors instead of the matrices, and generates the output directly from there.

% Parameters.
dim = 4;

% Test input.
x = [2, -2, 5, 4];
y = [1, -4, 6, -2];
z = rand(dim);
[X, Y] = meshgrid(x, y);
Z = griddata(x, y, z, X, Y);
[R, R_indx] = sort(Y);
[Q, Q_indx] = sort(X, 2);

% Initialize output.
S = zeros(dim);

% Provided solution using loop.
for i = 1:numel(z)
  for j = 1:numel(z)
    if (X(j) == Q(i) && Y(j) == R(i))
      S(i) = Z(j);
      break
    end
  end
end

% Output.
S

% Solution using array indexing; output.
S_array = reshape(((X(:) == Q(:).') & (Y(:) == R(:).')).' * Z(:), dim, dim)

% Solution using vector indexing; output.
[r, r_indx] = sort(y);
[q, q_indx] = sort(x);
[X, Y] = meshgrid(q, r);
Z = griddata(q, r, z, X, Y);
idx = (ones(dim, 1) * ((q_indx - 1) * dim) + r_indx.' * ones(1, dim));
S_vector = Z(idx)

Example output:

S = 
   0.371424   0.744220   0.777214   0.778058
   0.580353   0.686495   0.356647   0.577780
   0.436699   0.217288   0.883900   0.800133
   0.594355   0.405309   0.544806   0.085540

S_array =
   0.371424   0.744220   0.777214   0.778058
   0.580353   0.686495   0.356647   0.577780
   0.436699   0.217288   0.883900   0.800133
   0.594355   0.405309   0.544806   0.085540

S_vector =
   0.371424   0.744220   0.777214   0.778058
   0.580353   0.686495   0.356647   0.577780
   0.436699   0.217288   0.883900   0.800133
   0.594355   0.405309   0.544806   0.085540