I want to reduce a n-order ordinary differential equation into a first order system of equations. This is in preparation for numerical analysis. I use both Sympy and Sagemath for Computer Algebra, but I have not found any functions in them to do this type of reduction. I was not sure if anyone else could indicate whether this functionality exists within either Sympy or Sagemath.
An example of this would be reducing the equation:
x''' - 2x'' + x' = 0 to a system of first order equations:
[0 1 0 0 0 1 0 -1 2]