In this work, we propose an efficient and accurate numerical scheme for computing both the ground and first excited states of the nonlinear fractional Schrödinger eigenvalue problem. The method is based on a fractional gradient flow with discrete normalization, reformulated via the invariant energy quadratization (IEQ) approach to handle the nonlinear term in a linear way. Temporal discretization is performed using a second-order Crank-Nicolson scheme, while spatial discretization employs a Jacobi-Galerkin spectral method, yielding spectral convergence in space. We rigorously prove the discrete energy-decay property with respect to a modified energy functional. The resulting fully discrete scheme can compute different states simply by varying the initial condition, without altering the algorithmic framework. Numerical experiments confirm its high accuracy, efficiency, and robustness, and demonstrate its capability to capture intricate solution structures across a wide range of fractional orders and nonlinear interaction strengths.