Which analogy between polynomials and differential equations did Rota have in mind in his TEN LESSONS? Rota writes in his TEN LESSONS I WISH I HAD LEARNED BEFORE I STARTED TEACHING
DIFFERENTIAL EQUATIONS:

"It states that every differential polynomial in the two solutions of
  a second order linear differential equation which is independent of
  the choice of a basis of solutions equals a polynomial in the
  Wronskian and in the coefficients of the differential equation (this
  is the differential equations analogue of the fundamental theorem on
  symmetric functions, but keep it quiet)."

And:

"Worse, no one realizes that changes of variables are not just a
  trick; they are a coherent theory (it is the differential analogue of
  classical invariant theory, but let it pass)."

And:

"For second order linear differential equations, formulas for changes
  of dependent and independent variables are known, but such formulas
  are not to be found in any book written in this century, even though
  they are of the utmost usefulness. Liouville discovered a differential
  polynomial in the coefficients of a second order linear differential
  equation which he called the invariant. He proved that two linear
  second order differential equations can be transformed into each other
  by changes of variables if and only if they have the same invariant.
  This theorem is not to be found in any text. It was stated as an
  exercise in the first edition of my book, but my coauthor insisted
  that it be omitted from later editions."

Where can I learn more about this?
 A: This might be a partial answer, but I will add later if need be.
There is quite a bit of theory underlying a systematic change of variables.  The theory underlying it all is Noether's Theorem, which states for any symmetry there is a corresponding invariant.  This invariant then is a useful substitution, as it will effectively reduce the order of a differential equation by 1, and if the order is already 1, the substitution makes the equation separable.
The standard example of this is with first order homogenous equations.  In this case, we can see that the differential equation has symmetry under the transformation $x \rightarrow \lambda x$ and $y \rightarrow \lambda y$.  An invariant corresponding to the symmetry is some quantity that does not change under the transformation; thus, $y/x$ is an invariant.
There are many other examples of such a change of variable.  For example, take the differential equation
$$x^{3/2}y''+\sqrt{x}y'+\frac{y^2}{x}=1$$
In this case, you can check that that the differential equation is symmetric with respect to the transformation $x \rightarrow \lambda^2 x$, $y \rightarrow \lambda y$.  With this transformation, $y^2/x$ is an invariant of the transformation, and therefore the order of the equation can be reduced by 1 by substituting $u=y^2/x$.  (specifically, it changes the equation into a differential equation that contains $u''$ and $u'$, so a substitution can be made to reduce the order).
There are however many other symmetries, and corresponding invariants, possible however, beyond these "scaling" symmetries.  So the question then becomes how to find such substitutions.  The most complete and systematic method is via Lie Groups to calculate the symmetries.  Some references for doing so are here and here.  The most comprehensive resource on this question is here.  You will notice that these resources are generally graduate level however; I am not aware of systematic treatments of calculating these symmetries (beyond scaling symmetries) at a lower level.
A: I'm adding my own answer because I figured out a small piece of the puzzle. Perhaps other answers can build on this and maybe this makes it clearer what kind of thing I'm looking for.
When we have a polynomial equation like $x^3 + bx^2 + cx + d = 0$ we can express symmetric polynomials of the roots as polynomials in the coefficients. For instance, if $(x - r_1)(x - r_2)(x - r_3) = x^3 + bx^2 + cx + d$ then $r_1 r_2 r_3 + r_1 + r_2 + r_3 = -(a+d)$. This is cool because the roots might be algebraic numbers, but any symmetric polynomial in them can be calculated explicitly, and will be a rational number if the coefficients of the original polynomial are rational.
Using the Laplace transform we can transform a differential equation like $f''' + bf'' + cf' + df = 0$ to the previous polynomial equation. This will give us three linearly independent $f_1 = \exp(r_1 t), f_2 = \exp(r_2 t), f_3 = \exp(r_3 t)$, provided the roots are distinct.
Suppose that we have a term $r_1^n r_2^m r_3^k$ in the symmetric polynomial. We then do $(D^n f_1)(D^m f_2)(D^k f_3) = r_1^n r_2^m r_3^k \exp((r_1 + r_2 + r_3) t)$. Hence if we have a symmetric "differential polynomial" in the $f_1,f_2,f_2$ then that will be some symmetric polynomial in the roots times $\exp((r_1+r_2+r_3)t)$. Since both the symmetric polynomial in the roots and the $r_1+r_2+r_3$ can be found as polynomials in the coefficients $b,c,d$, we can explicitly calculate any symmetric differential polynomial of the three solutions.
Now, there are many questions:


*

*What if we picked another set of linearly independent solutions?

*What about repeated roots?

*What if the equation is not constant coefficients?


I think an anwser to these questions would explain what Rota meant by the first paragraph that I quoted. 
Lie groups might indeed have been what he meant by the second paragraph, but I'd like to more clearly understand how they are analogous to classical invariant theory.
And then there's the third paragraph...
I don't know if I can award a bounty multiple times, but if I can, then I will award an additional 50 point bounty for an explanation for each of those three paragraphs.
A: TL;DR: For the first two excerpts see
Kung, Joseph P. S.; Rota, Gian-Carlo
On the differential invariants of a linear ordinary differential equation.
Proc. Roy. Soc. Edinburgh Sect. A 89 (1981), no. 1-2, 111–123.
For the last excerpt see the exercises on page 57 of the 4th edition of the book Rota mentions.

The first excerpt refers to a theorem proved by Appell in 1881 (http://www.numdam.org/item/ASENS_1881_2_10__391_0/). According to the paper by Kung-Rota, the novelty of their paper is that they prove the theorem in a purely algebraic manner, whereas Appell uses analysis. I didn't look at Appell's paper, but I can verify that Kung-Rota uses no analysis. Their proof can be considered as a meditation on the universality of the determinant. Here is the abstract:

The second excerpt is related to a remark in the same paper on the elimination theory developed by the Ritt school. The analogy between classical invariant theory and linear differential equations seems to be about what is called "differential algebra", that is, algebra of $+,-,\times,\div,\;'$, the prime acting linearly and according to the product rule for derivatives.

For the third excerpt, consider the "spaces of $\mathcal{C}$oefficient functions"
$$\mathcal{C}:= C^1(\mathbb{R},\mathbb{R})\times C^0(\mathbb{R},\mathbb{R}), \quad\tilde{\mathcal{C}}:= C^0(\mathbb{R},\mathbb{R})\times C^1\left(\mathbb{R},\mathbb{R}_{>0}\right),$$
where $C^r$ means "still continuous after $r$ times differentiated". Each pair in either one of these spaces defines a linear second order ODE by
$$(p,q) \mapsto \quad u''+pu'+qu=0.$$
Observe that we can take both $p$ and $q$ to be only continuous for the classical existence and uniqueness theory. Differentiability is required to define the invariant $\mathbb{I}$ Rota refers to, provided that we don't use distributions. Similarly for $\mathbb{J}. $Define
$$\mathbb{I}:\mathcal{C}\to C^0(\mathbb{R},\mathbb{R}),\quad (p,q) \mapsto q-\dfrac{p^2}{4}-\dfrac{p'}{2}$$
and
$$\mathbb{J}:\tilde{\mathcal{C}}\to C^0(\mathbb{R},\mathbb{R}),\quad (p,q) \mapsto \dfrac{q'+2pq}{q^{3/2}}.$$
Birkhoff-Rota has these two as a collection of exercises:
Proposition (Change of dependent variable for $\mathcal{C}$): There is a function $\varphi\in C^2(\mathbb{R},\mathbb{R})$ such that $u$ solves the equation defined by $(p_1,q_1)$  iff $e^\varphi u$ solves the equation defined by $(p_2,q_2)$, i.e. the ODEs are $C^2$-equivalent, if and only if
$$\mathbb{I}(p_1,q_1) = \mathbb{I}(p_2,q_2).$$
Proposition (Change of independent variable for $\tilde{\mathcal{C}}$): There is a diffeomorphism $\varphi \in C^2(\mathbb{R},\mathbb{R})$ with $C^2$ inverse such that $u$ solves the equation defined by $(p,q)$ iff $u\circ \varphi^{-1}$ solves a constant coefficient second order linear ODE, i.e. the ODE is $C^2$-equivalent to a constant coefficient linear ODE, if and only if
$$\mathbb{J}(p,q)=\text{constant}.$$
The proofs of these are not complicated once reduction steps are known.
Let me point out that these are very explicit algebraic calculations, and there are no Baire-category type arguments involved. Finally observe that the last proposition gives a nice criterion for a second order linear ODE in a relatively large class to be explicitly solvable by second-year undergraduate methods.
