Suppose we want to convert, \begin{align*} y_1’’ &= t^2-y_1’-y_2^2 &\quad \quad y_1(0)=0\quad y_1’(0)=1\\ y_2’’ &= t-y_2’-y_1^3 &\quad \quad y_2(0)=1\quad y_2’(0)=0 \end{align*} into a first-order autonomous linear system of ODEs.
To convert to first-order and autonomous you can define, \begin{align*} x_1 &= y_1 \\ x_2 &= y_1’ \\ x_3 &= y_2 \\ x_4 &= y_2’ \\ x_5 &= t \end{align*} and we end up with \begin{align*} x_1’ &= x_2 \\ x_2’ &= x_5^2-x_2-x_3^2 \\ x_3’ &= x_4 \\ x_4’ &= x_5+x_4-x_1^3 \\ x_5’ &= 1 \end{align*} But it remains to convert it to a linear system. How does one accomplish that? Can one simply define $x_6=y_2^2$ and $x_7=y_1^2$ and $x_8=y_1^3$?
It is not so easy. The thing with monomials in a function is that it's differential: $$(y(t)^n)' = y'(t)\cdot ny(t)^{n-1}$$
This is product between function and it's derivative. It is non-linear.
I think you need to find some representation where function concatenation with polynomial is linear. This is possible with for example Carleman matrices.
So our function y is represented $\bf M_y$, squaring represented by $\bf M_2$, Differentiation $\bf M_D$
Now we can write $(y(t)^2)'$ as $$\bf M_D M_2M_y$$
But I am not sure that such $\bf M_D$ exists for these Carleman matrices.