It does not appear that the data you plot is the data you provide in the .mat file. If not, then how can we rationally help you?
The pictures that you show are of a function that is not apparently single valued, in the form z(x,y). The plots you show have a function that rolls over on top of itself. A single valued relation is crucial for griddata to work. (The same requirement would apply for my own gridfit. So don't try gridfit and then get upset when it fails to work well either.)
Note that the data you do provide appears to have a singularity at x=y=0. That will tend to make it difficult to produce a well behaved surface, due to the extremely high slopes.
Lets look at the data you have given us, plotted just as (x,y) pairs.
plot(x,y,'.')
MATLAB does not magically know how to connect these points. All you have given us is a list of those points as scattered pairs. So while your eye may be able to connect them into a surface, all that MATLAB can do is to form a delaunay triangulation in (x,y). It has no idea how the points are connected otherwise.
Here is a small portion of the triangulation that would have been formed:
These long, thin triangles are a problem in interpolation, as well as the choice of which triangles were formed. When you have a relationship with huge gradients, linear interpolation using long thin triangles is never going to work well. That, by the way, is what griddata uses here by default. Yes, it is true that the other methods of griddata will also work poorly on this data. Again, singularities produce problems. If singularities were easy to deal with, they would not have a special name for them, would they?
As well, the fact that x is scaled from 0 to 5e3 and y varies from 0 to 5e6 is also a problem, since it makes linear interpolation that is based on a delaunay triangulation difficult.
Worse, your data appears to be slightly concave along that top edge of the plot in (x,y). The result are some extremely long triangles along that edge. You can see them in the plot above. That will cause major problems.
One solution to the problem might be to rotate the entire data set in 3 dimensions, avoiding any issue of singularities, or single valued functions. That transformation might differ for each problem though, and you would need to understand how to do that.
A better solution would be to build the triangulation more intelligently. Then you could use tools to interpolate that functional form. But that would require you to learn to use such tools.
