In the paper, the multistep inverse finite element method (FEM) has been introduced to improve the accuracy of simulation in sheet metal stamping. Furthermore, the multistep inverse FEM can be used to obtain the strain/thickness distribution and shape of blank in the intermediate configurations. But there are three key problems, which are essential to implement multistep inverse FEM: the fist one is how to obtain the intermediate configurations of intermediate steps, the second one is how to find the corresponding coordinates in the sliding constraint surface, and the last one is how to update strain/stress distribution in the intermediate configurations in a fast and reliable way. Based on the known configurations of punch and die of the current step, the strategy of area minimization coupled with feasible sequential quadratic programming code is used to obtain initial intermediate configurations. An efficient walk-through point location algorithm with its complexity per point ( means the space dimension) is used to deal with contact searching problem and restrain the movement of corresponding nodes of intermediate configurations. In order to preserve the computational efficiency of inverse FEM, a pseudodeformation theory of plasticity based constitutive equation is proposed, which can well reflect the actual forming condition such as elastic/plastic deformation or loading/unloading condition. The above-mentioned improvements are implemented in our in-house inverse analysis software INVERSTAMP/MULTISTEP module. The presented algorithms are applied to a two-step cylinder cup deep-drawing product and three-step S-rail forming case. The numerical results compared with explicit dynamic solver LS-DYNA3D confirm its validity in formability prediction of intermediate shapes and final workpiece.