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