Hi Chris,
That lat and lon variables in your code span from [0 30] and [0 20] respectively, as opposed to the actual limits set on the axes by domain.lat and domain.lon. The scaling equations are probably not correct. Also, in surfacem, the fourth argument should be "ZData" and third corresponds to "CData". You can verify the above statements by using the following code
figure;
view(3)
ax1 = axesm('mapprojection','miller', ...
'frame','on');
setm(ax1, 'maplatlimit', [0 30], ...
'maplonlimit', [0 20], ...
'flatlimit', [0 30], ...
'flonlimit', [0 20]);
axis normal
tightmap
h1 = surfacem(lat, lon,z,zeros(size(z)));
zlim([-1 1])
Given x and y are symmetric around origin, you can use the following equations to calculate lon and lat.
lat = y.*(max(domain.lat) - min(domain.lat))/(max(y(:)) - min(y(:))) + mean(domain.lat);
lon = x.*(max(domain.lon) - min(domain.lon))/(max(x(:)) - min(x(:))) + mean(domain.lon);
