z = [-1:0.1:1]
theta = [0:pi/18:2*pi];
[Z, T] = meshgrid(z, theta)

R = sqrt((1-Z.^2)./(0.25*cos(T).^2 + sin(T).^2))
X = R.*cos(T)
Y = R.*sin(T)
mesh(X,Y,Z)
axis equal