I used to solve differential algebraic equations using Lagrange polynomials.

Essentially you convert the differential equations into an algebraic system by discretizing the solution. The method is called Orthogonal Collocation on Finite Elements (OCFE), and it was developed by chemical engineers.

The Lagrange polynomials were calculated at special knots that corresponded to Radau interior points, which work great for stiff systems.

It’s great for solving differential algebraic equations through purely sparse matrix operations, no explicit integration like Runge Kutta. (Well, it’s implicit Runge Kutta).