Simulation of stiff multibody systems requires Jacobian matrices (JM) of equations of motion. Evaluation of the JM by finite differences is a very CPU time-consuming operation. Use of approximate JM taking into account stiff forces reduces considerably the computational efforts. Analytic expressions for the corresponding matrices are obtained. Block-diagonal approximations of the JM are introduced to apply implicit solvers to the scheme of the articulated body algorithm as well as to simulate system of thousands of bodies undergoing contact interactions. Models of a fright coach and a ballast system illustrate implementation of the developed approaches.