Without looking deeply into the method used to interpolate, I can't absolutely know what their code does. However, I do deeply understand interpolation in 1 and more dimensions, and what is commonly done. (In a previous incarnation, in a galaxy far, far away, I spent a fair amount of time learning about interpolation artifacts in multiple dimensions, why they arise, etc.)
First, let me point out that linear interpolation on a problem with two independent variables (so z(x,y)) requires three points. That is, THREE points determine a plane. And a linear interpolation implicitly assumes a plane. So what happens? What is done when you have more than 3 points in a rectangular lattice? After all, a rectangle has FOUR points, not three. The difference is huge. Well, it can be significant.
Consider a simple example in two dimensions, on the unit square [0,1]X[0,1]. I'll pick some values for the 4 corners of that square, where there is an intentional nonlinearity. No single plane will pass through all 4 points.
Fij = [0 3;1 1.5]
First, I'll plot that array using what people often call a bilinear interpolant.
F = @(x,y) Fij(1,1)*(1-x).*(1-y) + Fij(1,2)*(1-x).*y + Fij(2,1)*x.*(1-y) + Fij(2,2)*x.*y;
fsurf(F,[0 1 0 1])
And what you see is a curved surface. Again, nothing you can do will make that curved thing into a plane. And while you talked about doing "linear" interpolation, I think tyou do not really understand linear interpolation, and what happens in more than one dimension. What the function F implements is what is sometimes called bilinear interpolation. It truly is linear if you look along any line parallel to either axis, but that is as much as you can do.
Instead, the common approach is to dissect each square in a rectangular grid into a pair of triangles. Perhaps something like this:
trimesh([1 2 4;1 3 4],[0;0;1;1],[0;1;0;1],[0;3;1;1.5])
The result is two planes.
trimesh([1 2 3;2 3 4],[0;0;1;1],[0;1;0;1],[0;3;1;1.5])
In all cases, the 4 corners of that square have identical function values. But we get three different surfaces. In cases 2 and 3, the only difference is how I oriented the triangular dissection of the unit square.
Now, what happens when we plot a surface on a grid? MATLAB will choose a specific dissection, breaking EACH rect in that grid into two triangles. I don't know which way the dissection is done internally of course. It might be this one:
Or it might be one where the triangles are oriented 90 degrees to that, thus...
What matters is the choice made, and I cannot see into their internal code. But the choice made will be consistent internally. It will be one of those two alternatives.
But what happens when you flip the data along one axis? Implicitly the triangulation is flipped compared to your data. They are still using the same internal dissection. But your data is flipped.
And that means the interpolation used to interpolate the color shading changes, in a subtle way. As you can see, those 2nd and 3rd surfaces are actually pretty significantly different. Enough so that when you look at the resulting coloration as generated, we will see the differences, because those differences, subtle as they may seem, will impact where the color quantization changes.
Now, let me look at your specific data.
data = [0.489352000000000 0.484068000000000 0.511895000000000 0.529349000000000 0.550050000000000 0.589397000000000 0.607263000000000 0.641485000000000 0.662575000000000 0.687064000000000 0.713925000000000
0.507529000000000 0.527680000000000 0.538550000000000 0.557169000000000 0.577756000000000 0.614422000000000 0.657682000000000 0.676427000000000 0.684739000000000 0.714883000000000 0.738783000000000
0.538647000000000 0.548932000000000 0.566195000000000 0.580673000000000 0.618064000000000 0.654818000000000 0.689324000000000 0.702041000000000 0.719743000000000 0.733533000000000 0.764315000000000
0.556159000000000 0.573411000000000 0.596687000000000 0.623561000000000 0.671057000000000 0.689452000000000 0.714168000000000 0.730792000000000 0.748139000000000 0.774988000000000 0.788946000000000
0.590044000000000 0.600503000000000 0.634420000000000 0.648230000000000 0.685438000000000 0.721476000000000 0.746399000000000 0.759050000000000 0.782877000000000 0.790311000000000 0.816816000000000
0.620506000000000 0.638116000000000 0.665363000000000 0.676466000000000 0.716868000000000 0.749803000000000 0.754542000000000 0.779754000000000 0.801076000000000 0.827632000000000 0.854378000000000
0.663773000000000 0.672726000000000 0.696712000000000 0.710226000000000 0.735066000000000 0.747308000000000 0.775699000000000 0.807429000000000 0.828759000000000 0.855226000000000 0.875827000000000
0.695315000000000 0.709644000000000 0.726708000000000 0.740776000000000 0.768420000000000 0.799622000000000 0.796054000000000 0.825712000000000 0.846637000000000 0.876592000000000 0.897063000000000
0.719426000000000 0.727156000000000 0.748276000000000 0.762108000000000 0.796051000000000 0.804258000000000 0.813496000000000 0.820555000000000 0.868153000000000 0.888011000000000 0.908845900000000
0.743718000000000 0.754901000000000 0.765970000000000 0.776975000000000 0.811003000000000 0.825777000000000 0.828351000000000 0.841359000000000 0.869123000000000 0.903499400000000 0.914743600000000
0.765627000000000 0.770041000000000 0.787120000000000 0.801125000000000 0.823051000000000 0.837764000000000 0.845964000000000 0.862046000000000 0.883647000000000 0.895037000000000 0.919736800000000];
x = 0:10;
y = 0:10;
[xgrid, ygrid] = meshgrid(x,y);
A big part of the problem, what probably exacerbates the issure, is you are trying to use a highly quantized colormap, thus jet(32). If you had used a finer colormap, perhaps jet(256), the differences would be far less obvious.
figure
colormap(jet(256))
surface(xgrid,ygrid, data, EdgeAlpha=0.2, FaceColor='interp');
