In this study, we present a Finite Difference Method (FDM)-based stress integration algorithm for Crystal Plasticity Finite Element Method (CPFEM). It addresses the complexity of computing the first derivative of resolved shear stress in the Euler backward stress integration algorithm with Newton-Raphson method. The proposed FDM-based model was verified by evaluating its accuracy, convergence and computational efficiency through single-element simulations. The developed FDM-based model can be easily applied to various constitutive models for CPFEM, overcoming the problem of deriving complex derivative regardless of constitutive models. Additionally, the proposed FDM-based model was validated with the reduced texture approach using AA 2090-T3. Specific parameters including crystallographic orientations were calibrated and the plastic anisotropy was successfully described. In addition, the earing profiles were compared using various stress integration methods. As a result, the proposed FDM-based model can be used as an alternative to the Euler backward method using analytic derivatives with the compatible accuracy, convergence, computational efficiency along with easy implementation within the CPFEM framework.