This study deals with stress update algorithm based on finite difference method. The developed algorithm based on multi-stage Euler backward method where the first and second derivatives of yield function are approximated by central difference method. With the developed algorithm, it is possible to conduct elastic-plastic finite element simulation without analytical first and second derivatives of yield function, which has been the biggest obstacle when advanced constitutive models such as Yld2000-2d, Yld2004, and homogeneous anisotropic hardening (HAH) model are implemented for finite element modeling. For the verification purpose of the developed algorithm, single element simulation for r-value and stress directionalities prediction was conducted with Hill48 and Yld2000-2d anisotropic yield function under associated and non-associated flow rules. The simulation results were compared with the theoretical predictions and the results obtained from analytical Euler backward and Euler forward methods. To check the availability on distortional plasticity such as HAH model, single element tension followed by compression and tension (RD) followed by tension (TD) simulation were carried out and compared with the reference data. Finally, simulation for earing prediction with the developed algorithm was performed with various advanced constitutive models: Hill48 and Yld2000-2d under both associated and non-associated flow rules. The simulation results show a strong availability of the developed algorithm as an alternative to classical Euler backward method.