This can also be solved using the patch() function to plot the parametric equation of a circle:
% Create and plot sphere grid
[X,Y,Z] = meshgrid(-2.5:.5:2.5,-2.5:.5:2.5,0:.5:5);
scatter3(X(:),Y(:),Z(:),'filled')
% Create and plot circle
C = [0,0,0]; % Center of circle
R = 1 ; % Radius of circle
theta = 0:0.01:2*pi;
x = C(1)+R*cos(theta);
y = C(2)+R*sin(theta);
z = C(3)+zeros(size(x));
hold on
patch(x,y,z,'k');
hold off
axis equal
