This is an interpolation artifact. The cause is a pretty simple one. Fixing it is slightly less trivial, because you need to understand what happened, AND you need to understand interpolation sufficiently to avoid the artifact.
First, 3 points determine a plane. Second, in order to perform truly linear interpolation for a problem z(x,y), you need a planar surface. And what is rule # 1? 3 points determine a plane.
However, the arrays X,Y,Z are rectangular things. In a plot of 3x3 arrays like this, there are 4 rectangular cells. I'll pick just one of them for a closeup.
format short g
Recognize that this little 4 cornered square is composed of 4 points. FOUR. Not THREE.
So what should MATLAB do? You cannot linearly interpolate across that domain. You might consider what is sometimes called bilinear interpolation, but that too introduces artifacts. What is typically done is to break the region down into TWO triangles.
Then interpolate across each triangle, now using LINEAR interpolation. 3 points allow linear interpolation. But what happens?
Remember that your array is REALLY coarse. We see an interpolation artifact, because we can effectively see the shared edge of the pair of triangles the square was dissected into.
Can that be improved? Perhaps. The improvement would require that you more intelligently interpolate the data. By a more fine interpolation of the array, and doing so smoothly, we reduce any visible interpolation artifacts. By the way, there can still be artifacts generated, but they will be REALLY TINY now, effectively so small as to be invisible. Don't allow surf to make the decision about how to interpolate that very coarse array. If you do, then you get these highly visible interpolation artifacts.
[Xint,Yint] = meshgrid(linspace(-0.0015,0.0015,101),linspace(-0.0005,0.0005,101));
Zint = interp2(X,Y,Z,Xint,Yint,'makima');
H = surf(Xint,Yint,Zint,'FaceColor','interp');
H.EdgeColor = 'none';
asp_ratio = daspect; daspect([1 1 asp_ratio(3)]); colormap(jet(10));
H2 = mesh(X,Y,Z);
H2.FaceColor = 'none';
H2.EdgeColor = 'k';
I've added the mesh plot in black so you can see the original array "edges".
At the end, the only artifacts I see are those imposed by a very short colormap, thus jet(10). jet(256) results in a fairly reasonable rainbow colored surface.
As an alternative, I could have used what interp2 calls 'linear' interpolation.
Zint = interp2(X,Y,Z,Xint,Yint,'linear');
In fact, that is something often called bilinear interpolation, or perhaps tensor product linear interpolation.
This may be close to what you had drawn as the desired surface.