We present a three-step magnetic inversion procedure in which invariant quantities with respect to source parameters are inverted sequentially to give (1) shape cross section, (2) magnetization intensity, and (3) magnetization direction for a 2D (elongated) magnetic source. The quantity first inverted (called here the shape function) is obtained from the ratio of the gradient intensity of the total-field anomaly to the intensity of the anomalous vector field. For homogenous sources, the shape function is invariant with source magnetization and allows reconstruction of the source geometry by attributing an arbitrary magnetization to trial solutions. Once determined, the source shape is fixed and magnetization intensity is estimated by fitting the total gradient of the total-field anomaly (equivalent to the amplitude of the analytic signal of magnetic anomaly). Finally, the source shape and magnetization intensity are fixed and the magnetization direction is determined by fitting the magnetic anomaly. As suggested by numerical modeling and real data application, stepped inversion allows checking whether causative sources are homogeneous. This is possible because the shape function from inhomogeneous sources can be fitted by homogeneous models, but a model obtained in this way fits neither the total gradient of the magnetic anomaly nor the magnetic anomaly itself. Such a criterion seems effective in recognizing strongly inhomogeneous sources. Stepped inversion is tested with numerical experiments, and is used to model a magnetic anomaly from intrusive basic rocks from the Paraná Basin, Brazil.