It seeks a highly accurate sparse identification method for the nonlinear dynamical systems, combinating with traditional numerical analysis techniques. Firstly, a suitable basis function library must be created in order to approximate the potential nonlinear dynamical systems. Then, the approximated nonlinear dynamical systems are discretized using the linear multistep method. Next, when the state data contains noise, the generalised least squares method is used to calculate the approximate covariance matrix of the noise residual term and use this matrix to weight the optimisation problem obtained from the preceding process, thereby reducing the influence of noise on the model identification results. Finally, the subspace pursuit algorithm selects the set of features with the smallest coefficient error from the data to serve as the basis function library for the next iteration, and after the iteration is completed, the coefficient values of the retained features are computed using the least squares method. The proposed linear multistep subspace pursuit methods for identifying nonlinear dynamic systems possess high accuracy and robustness. Numerical results are presented to demonstrate the effectiveness of the proposed methods.