We have developed a gravity inversion method to estimate a discontinuous basement relief based on an extended version of Bott’s method that allows variable density contrasts between sediments and basement, optimizes the modulus of the solution correction at each iteration, and provides for solution stabilization. Initially, we obtain a linear approximation stabilized by the total variation functional that correctly maps the horizontal positions of the existing high-angle faults but produced poor estimates of the basin depths. Subsequent iterations update the depth estimates toward the correct values, at the same time preserving the correct fault horizontal positions. Additionally, we stabilize each solution correction by the smoothness constraint without inverting any matrix. The method was substantially more efficient than the nonlinear method, which solves a system of linear equations by the conjugate gradient method at each iteration. For 3000 parameters, it is almost four times faster than the conjugate gradient, and this ratio increased with the number of parameters.