Try using meshgrid() to get all possible combinations of x and y in a range you specify:
% Declare some constants.
fontSize = 20;
a=1;
b=2;
c=3;
numElements = 30;
constant = 1; % or 0, whichever you want.
% Define range of x and y axes.
x = linspace(-2*a, 2*a, numElements);
y = linspace(-2*b, 2*b, numElements);
% Get every possible combination of x and y.
[X, Y] = meshgrid(x, y);
% Construct function.
Z1 = c * sqrt(constant + X.^2 / a^2 + Y.^2 / b^2);
Z1 = reshape(Z1, [], length(y));
surf(Z1);
Z2 = -c * sqrt(constant + X.^2 / a^2 + Y.^2 / b^2);
Z2 = reshape(Z2, [], length(y));
hold on;
surf(Z2);
% Fancy up the plot
grid on;
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
zlabel('Z', 'FontSize', fontSize);
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
