In this paper, we propose a method based on a two-stage algorithm to simultaneously identify the coefficients and fractional differentiation orders of fractional order systems (FOSs) with commensurate order. The proposed method adopts the fractional integral operational matrix of block pulse functions (BPFs) to convert the FOS to a linear parameter regression equation. Then, a two-stage algorithm is developed to identify the coefficients and orders. First, with the orders fixed, the coefficients are identified using the instrumental variable-based recursive least square algorithm. Then, with the identified coefficients fixed, the orders are estimated using the Gauss–Newton iterative algorithm. The above process iterates until the stop criterion is met. Two identification examples are given to verify the effectiveness of the proposed method.