In this work, a framework for the automatic symbolic index reduction and numerical integration of generic differential-algebraic equation systems is presented. The proposed approach does not exploit any a priori knowledge of the specific differential-algebraic system of equation formulation. The index reduction is performed only with the aim of linear algebra techniques. Hierarchical representation of expressions is conveniently used to limit expression swell and ensure the numerical stability of the solution. The effectiveness of the algorithm is validated through symbolic and numerical experiments on a multi-body model.