In this paper we investigate the behavior of numerical ODE methods for the solution of systems of differential equations coupled with algebraic constraints. Systems of this form arise frequently in the modelling of problems from physics and engineering; we study some particular examples from electrical networks, fluid dynamics and constrained mechanical systems. We show that backward differentiation formulas converge with the expected order of accuracy for these systems.