Considering the flexibility and variety of the distributed generation system structure， the numerical differentiation method is adopted to get the Jacobian matrix for implicit integration algorithm. Then a numerical integration scheme suitable for stability simulation of the distributed generation system is constructed subsequently. Through this algorithm， Jacobian matrix acquisition process can be simplified and the procedure exhibits parallel characteristics. Moreover， the proposed method is especially suitable for stability simulation involving variable-structure models. Simulation results are given to validate the accuracy and efficiency of presented algorithm.