M = 1;
Q = 1;
horizon = M + sqrt(M*M - Q*Q);
dr = 0.2;
r = [horizon+dr:dr:10*horizon];
dzdr = sqrt((2*M*r - Q*Q) ./ (r.*r - 2*M*r + Q*Q));
z = cumsum(dzdr * dr);

theta = [0:pi/16:2*pi];

rr  = transpose(r) * ones(size(theta));
TT  = transpose(ones(size(r))) * theta;
zz  = transpose(z) * ones(size(theta));

rrr = transpose(r([length(r):-1:1])) * ones(size(theta));
rzz = transpose(z([length(z):-1:1])) * ones(size(theta));

rr = [ rrr; rr];
TT = [  TT; TT];
zz = [-rzz; zz];

xx = rr .* cos(TT);
yy = rr .* sin(TT);

surf(xx,yy,zz);
