We present a mathematical and numerical investigation to the shrinking-dimer saddle dynamics for finding any-index saddle points in the solution landscape. Due to the dimer approximation of Hessian in saddle dynamics, the local Lipschitz assumptions and the strong nonlinearity for the saddle dynamics, it remains challenges for delicate analysis, such as the boundedness of the solutions and the dimer error. We address these issues to bound the solutions under proper relaxation parameters, based on which we prove the error estimates for numerical discretization to the shrinking-dimer saddle dynamics by matching the dimer length and the time step size. Furthermore, the Richardson extrapolation is employed to obtain a high-order approximation. The inherent reason of requiring the matching of the dimer length and the time step size lies in that the former serves a different mesh size from the later, and thus the proposed numerical method is close to a fully-discrete numerical scheme of some space-time PDE model with the Hessian in the saddle dynamics and its dimer approximation serving as a "spatial operator" and its discretization, respectively, which in turn indicates the PDE nature of the saddle dynamics.