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 $y'(t) = Ay$ with

$$ A = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 2 \end{pmatrix}. $$