Factoring $n-$th order linear ode I'm trying to solve problem 3.36 of Gerald Teschl's ODE and Dynamical Systems, which asks to show that any linear $n-$th order equation can be factorized into first order equations:

Let $\phi_1,...,\phi_n$ be linearly independent solutions to the $n-$th order equation $L_n(f)=0$. Set
$$L_1(f):=\frac{W(\phi_1,f)}{\phi_1}=f'-\frac{\phi_1'}{\phi_1}f$$
and define $\psi_j:=L_1(\phi_j)$. Show that $\psi_2,...,\psi_n$ are linearly independent and
$$L_n(f)=L_{n-1}(L_1(f)),\qquad L_{n-1}(f)=\frac{W(\psi_2,...,\psi_n,f)}{W(\psi_2,...,\psi_n)}$$

where $W$ is the wronskian of functions.
I'm at a loss for how to proceed.
For linearly independence, suppose $$0=c_2\psi_2+...+c_n\psi_n=\phi_1(c_2W(\phi_1,\phi_2)+...+c_2W(\phi_1,\phi_n))$$
I would like to say that since $\phi_1,...,\phi_n$ are linearly independent, the pairwise Wronskians must also be, but I can't prove this (or don't see why its evident).
And similarly, I'm not sure how to proceed with the next parts. Any help is appreciated.
 A: For any $n$ functions $y_1,...,y_n$, the expression
$$
0=W(y_1,...,y_n,y)
$$
is a linear differential equation for $y$ of order $n$. Its leading coefficient is $W(y_1,...,y_n)$, so that the quotient in $L_{n-1}(f)=0$ produces a normalized linear differential equation that has $\psi_2,...,\psi_n$ as solutions.
As $L_1(\phi_1)=0$ and $L_2(\phi_k)=\psi_k$, $k=2,...,n$, the equation $L_{n-1}(L_1(f))=0$ has all of $\phi_1,\phi_2,...,\phi_n$ as solutions. As these functions are independent and the operator combination still has leading coefficient $1$, this completely determines the operator. Thus it is identical to any other operator with that order and solution basis that has leading coefficient $1$. This is most likely somewhere stated or strongly implied for $L_n$.

This whole construction formalizes the reduction-of-order method. Take as example some third-order equation $L_3[y]=y'''+c_2y''+c_1y'+c_0y=0$. If somehow a solution $\phi_1$ is found, then this method tries to find other solutions in the form $y=u\phi_1$. Inserting per product rule gives
\begin{align}
0&=[u'''\phi_1+3u''\phi_1'+3u'\phi_1''+u\phi_1''']
+c_2[u''\phi_1+2u'\phi_1'+u\phi_1'']+c_1[u'\phi_1+u\phi_1']
+c_0u\phi_1
\\
&=u'''\phi_1+[3\phi_1'+c_2\phi_1]u''+[3\phi_1''+2c_2\phi_1'+c_1\phi_1]u'
\end{align}
The terms containing $u$ cancel as $\phi_1$ is a solution. To get leading coefficient 1 one could divide by $\phi_1$. But that would give an extra coefficient in the operator composition identity. Instead absorb the leading coefficient into a new variable, set $\psi=u'\phi_1$ which also begets the formula of $\psi=L_1[y]$. Then the equation transforms to
\begin{align}
0&=\psi''+[\phi_1'+c_2\phi_1]u''+[2\phi_1''+2c_2\phi_1'+c_1\phi_1]u'
\\
&=\psi''+\left[\frac{\phi_1'}{\phi_1}+c_2\right]\psi'
+\left[2\phi_1''-\frac{\phi_1'^2}{\phi_1}+c_2\phi_1'+c_1\phi_1\right]u'
\\
&=\psi''+\left[\frac{\phi_1'}{\phi_1}+c_2\right]\psi'
+\left[2\frac{\phi_1''}{\phi_1}-\frac{\phi_1'^2}{\phi_1^2}+c_2\frac{\phi_1'}{\phi_1}+c_1\right]\psi
\\
&=L_2[\psi]
\end{align}
This shows that there is no easy structural relation between the coefficients of the original differential operator $L_3$ and the reduced $L_2$
