Like everything, there are probably better ways, especially to get the edges closed. The simple way is to just plot the surfaces and truncate them. Since the mesh is rectangular, so are the truncated edges.
% parameters
r = 2;
n = 1000; % number of points per axis
% only need one cylinder
[x1 y1] = meshgrid(linspace(-r,r,n));
z1 = sqrt(r^2 - x1.^2);
z2 = -sqrt(r^2 - x1.^2);
% truncate it
m1 = x1.^2 + y1.^2 > (r+r/n)^2;
z1(m1) = NaN;
z2(m1) = NaN;
% plot it and then plot a flipped copy
opts = {'specularstrength',0.5,'ambientstrength',0.4,'diffusestrength',0.9};
surf(x1,y1,z1,opts{:}); hold on;
surf(x1,y1,z2,opts{:})
surf(x1,z1,y1,opts{:})
surf(x1,z2,y1,opts{:})
axis equal
shading flat
colormap(ccmap)
lightangle(240,60)
lighting gouraud
