Stress integration algorithm based on finite difference method (FDM) was proposed to effectively deal with both first and second derivatives of yield and potential functions which are the lengthiest component in stress integration procedure. With the proposed numerical algorithm, both first and second derivatives of yield function are approximated by central difference method, so that finite element modeling using advanced constitutive model could be easily performed no matter how complicated its derivatives are. For the verification purpose, the algorithm was applied for advanced constitutive models: Plastic anisotropy model under associated (AFR) and non-associated flow rule (non-AFR), the homogeneous anisotropic hardening (HAH) model under associated (AFR) and non-associated flow rule (non-AFR). The proposed algorithm was verified with single element loading-unloading and cup-drawing simulations. The Euler backward method based on both the proposed numerical algorithm and analytical derivatives were employed for verification purpose. The accuracy and time efficiency of the proposed numerical algorithm were evaluated by comparing the simulation results from analytical derivatives. In addition, the applicability of the proposed numerical algorithm for the HAH models was estimated with single element loading-unloading, loading-reloading, and deep-drawing/springback simulations. Non-associated flow plasticity for the HAH model is newly proposed to improve numerical efficiency with finite difference method by keeping the same level of accuracy as associated flow rule plasticity. All the simulation results proved that the proposed numerical algorithm can be widely used for the implementation of advanced constitutive models. (C) 2018 Elsevier B.V. All rights reserved.