Actually, scatteredInterpolant function is easy to implement what you think.
function plot3D(x,y,z,v,m)
% scatter3(x,y,z);
Tmin = min(v); Tmax = max(v);
radius = mean(sqrt(x.^2+y.^2)); % use mean radius as radius
theta = linspace(0, 2*pi,m+1); % divide theta to m segments
minZ = min(z); maxZ = max(z);
Z = linspace(minZ, maxZ, m+1); % divide z to m segments
[Z, T] = meshgrid(Z, theta); % create meshgrid with theta and z
X = radius.*cos(T); % calculate X and Y for mesh
Y = radius.*sin(T);
F = scatteredInterpolant(x,y,z,v,'linear'); % create interp handle
V = zeros(size(X));
V(:) = F(X(:),Y(:),Z(:)); % interp the value
Vave = 0.5*(V(end,:)+V(1,:)); % fix the minor flaw (average the beginnig and the end of corresponding theta)
V(end,:) = Vave; % fix the minor flaw
V(1,:) = Vave; % fix the minor flaw
surf(X,Y,Z,V,'edgecolor','none')
caxis([Tmin Tmax]);colormap('jet');shading interp;colorbar;
Now use this function to produce figure
load('xyzvm_cylinder.mat')
plot3D(x,y,z,v,40)
The result shows as following
