Nevermind, I've figured out what I need. Something like this seems to work:
npts=100;
x = linspace(-1,1,npts);
y = linspace(-1,1,npts);
[X,Y]=meshgrid(x,y);
Z=X;
for r = 1:npts
for c = 1:npts
if (X(r,c)^2+Y(r,c)^2)>1
Z(r,c) = NaN;
else
Z(r,c)=sqrt(1-X(r,c)^2-Y(r,c)^2);
end
end
end
surf(X,Y,Z);shading interp
alpha(1)
[u v] = gradient(Z);
h = streamslice(X,Y,-u,-v);
set(h,'color','k')
for i=1:length(h);
zi = interp2(X,Y,Z,get(h(i),'xdata'),get(h(i),'ydata'));
set(h(i),'zdata',zi);
end
axis tight
%hold on;
%surf(X,Y,-Z);shading interp
%hold off;
daspect([1,1,1])
axis tight;