Actually, based on the sample picture you provided, you have a total of 17x21 grid points, and total number of cells in the grid is 16x20. You can use 'patch' function to obtain visualization similar to your sample image. Here is an example:
% Generate a M-by-N grid
[M,N]=deal(16,20);
x=linspace(-3,3,N+1);
y=linspace(-3,3,M+1);
[X,Y]=meshgrid(x,y);
V=[X(:) Y(:)]; % xy-coordinates of the grid-points
% Define connectivity of the inidividual cells making-up the grid
Np=numel(X); % total number of grid points
id=reshape(1:Np,[M+1 N+1]); % grid-point indices
i1=id(1:M,1:N);
i2=id(2:(M+1),1:N);
i3=id(2:(M+1),2:(N+1));
i4=id(1:M,2:(N+1));
F=[i1(:) i2(:) i3(:) i4(:)];
% Suppose these are the values of a scalar field measured at the grid points
Q=sinc(sqrt((X/2).^2+Y.^2))+sinc(sqrt((X).^2+(Y/3).^2));
Q=Q(:);
% Simulated grid currently lies in the xy-plane. Lets rotate it to make it
% look more like the grid in your image.
[~,Rz,~,Rx]=xyzRmat(pi*[45 0 30]/180);
V(:,3)=0;
V2=(Rx*Rz*V')';
% Here is one visualization appoach; values at the grid points will be
% interpolated to the interior of the cells
figure('color','w')
subplot(1,2,1)
patch('Faces',F,'Vertices',V2,'EdgeColor','k','FaceVertexCData',Q,'FaceColor','interp');
set(gca,'box','on')
axis tight
view([30 30])
camlight('headlight'), lighting phong
colorbar('Location','NorthOutside')
% Here is another visualization appoach; each grid cell will have a
% constrant value, corresponding to the average of its vertices
subplot(1,2,2)
patch('Faces',F,'Vertices',V2,'EdgeColor','k','FaceVertexCData',Q,'FaceColor','flat');
set(gca,'box','on')
axis tight
view([30 30])
camlight('headlight'), lighting phong
colorbar('Location','NorthOutside')
'xyzRmat' function used in this example is included as an attachment.
