How do I solve this ODE: $\frac{d^N y}{dx^N}=x^L y^M$ I have found a way to solve some unspecified equations of this type using Lie Theory:
$$\frac{d^N y}{dx^N}=x^L y^M$$ 
The answer is
$y=\bigg[\beta (\beta -1)(\beta-2)...(\beta-N+1)\bigg]^{\frac{1}{M-1}}x^{\beta}$ with $\beta = \frac{L+N}{1-M}$ 
I'll post the step-by-step solution in a few days, but I was wondering: aside from group theory, does anybody know another way to solve this problem?

I promised I would post the solution, so here it is.
In 1886 Sophus Lie stated that if a differential equation is invariant under the transformation of an infinite continuous group (Lie group), it can be rewritten as a function of the stabilizers of that group.  Lie also discovered that the relationships between the stabilizers could be found by taking their derivatives, and that these relationships, expressed in canonical form, could be used to solve the DEQ.  Dresner, Olver and others have noted that when the general solution cannot be derived, a special solution may often be found algebraically.  It was this that prompted me to propose and solve the above equation, and having solved it I began to wonder whether or not anybody else had heard of the technique. (I've used it before, but never on anything so generalized as this DEQ.)
Allow primed values to be the transformed variables, such that $x'=\lambda x$ and $y'=\lambda^\beta y$.  The unit conversion for this group occurs when $\lambda =1$. (Values of x' form a subgroup of y', and proving that all values of y' satisfy group properties is a trivial exercise.) Then $dx'=\lambda dx$, $dy'=\lambda^\beta dy$.  The transformed DEQ is 
$$\frac{d^N (\lambda^\beta y)}{(\lambda dx)^N}=(\lambda x)^L (\lambda^\beta y)^M$$
For invariance, $\lambda^{\beta -N}=\lambda^{L+\beta M}$ or $\beta = \frac{L+N}{1-M}$.  Once $\beta$ is found, invariance is assured.   
To find the stabilizers, first find the coefficients of the infinitesimal transformations and then use the method of characteristics.  For convenience, use Newton's notation such that $\dot{y}=\frac{dy}{dx}$, $\ddot{y}=\frac{d^2 y}{dx^2}$, etc.  Then $\dot{y}'=\lambda^{\beta -1}\dot{y}$, $\ddot{y}'=\lambda^{\beta -2} \ddot{y}$, $\dddot{y}'=\lambda^{\beta -3}\dddot{y}$ ....
$\frac{\partial x'}{\partial \lambda}\big|_{\lambda_o
=1}=x$, $\frac{\partial y'}{\partial \lambda}\big|_{\lambda_o
=1}=\beta y$, $\frac{\partial \dot{y}'}{\partial \lambda}\big|_{\lambda_o
=1}=(\beta-1) \dot{y}$, $\frac{\partial \ddot{y}'}{\partial \lambda}\big|_{\lambda_o
=1}=(\beta-2) \ddot{y}$, $\frac{\partial \dddot{y}'}{\partial \lambda}\big|_{\lambda_o
=1}=(\beta-3) \dddot{y}$, ..... 
$$d\lambda=\frac{dx}{x}=\frac{dy}{\beta y}=\frac{d\dot{y}}{(\beta -1) \dot{y}}=\frac{d\ddot{y}}{(\beta -2) \ddot{y}}=\frac{d\dddot{y}}{(\beta -3) \dddot{y}}=....$$ Integrating $\frac{dx}{x}=\frac{dy}{\beta y}$ produces $S_o =\frac{y}{x^\beta}$ where $S_o$ is a constant of integration, and constants are stabilizers for the group of polynomials. This is easily shown to be true by applying the group transformation to the stabilizer: it does not change.
$$S_o'=\frac{y'}{(x')^\beta}=\frac{\lambda^\beta y}{\lambda^\beta x^\beta}=S_o$$
Likewise, $\frac{dx}{x}=\frac{d\dot{y}}{(\beta -1) \dot{y}} \rightarrow S_1 =\frac{\dot{y}}{x^{\beta -1}}$, another stabilizer.  The next one is $\frac{dx}{x}=\frac{d\ddot{y}}{(\beta -2) \ddot{y}} \rightarrow S_2 =\frac{\ddot{y}}{x^{\beta -2}}$.  According to Lie there are an infinite number of stabilizers following the same pattern, so we can assume that $$S_{N-1}=\frac{\overset{(N-1)\bullet}{y}}{x^{\beta -N+1}}$$ $$S_N=\frac{\overset{N \bullet}{y}}{x^{\beta -N}}$$ 
Unless you are in error, at this point the DEQ may be rewritten in stabilizer form.  This is accomplished by dividing the left side by $x^{\beta -N}$ and dividing the right side by $x^{L+\beta M}$ (remembering that $\beta -N = L+ \beta M)$.  $$\frac{\overset{N \bullet}{y}}{x^{\beta -N}}=\frac{x^L y^M}{x^L y^{\beta M}}=\big(\frac{y}{x^\beta}\big)^M$$ $$S_N=S_o^M$$
Lie asserts that by taking the derivatives of these stabilizers we can find the canonical relationships between them.  With a little thought it is easy to see that $$x\frac{dS_o}{dx}=S_1-\beta S_o$$ $$x\frac{dS_1}{dx}=S_2 -(\beta -1)S_1$$ $$....$$  $$x\frac{dS_{N-1}}{dx}=S_N - (\beta -N+1)S_{N-1}$$
Within the direction field of the DEQ, singularities, saddle points and the separatrices that connect them are found in places where the above derivatives are all set equal to zero.  It is here that the special solutions are found.  Thus, $$S_1=\beta S_o$$ $$S_2=(\beta -1)S_1=\beta (\beta -1)S_o$$ $$....$$ $$S_N=\beta (\beta -1)(\beta -2)...(\beta -N+1)S_o$$
Substitute this last equation into the transformed DEQ and with a little algebraic effort you end up at the given solution.  
I really enjoy using this method because it exploits a group property of polynomials themselves, of which DEQ's are a subgroup.  It also turns the traditional numerical "shooting method" into a back-of-the-envelope calculation. 
Sophus Lie was a genius.  I highly recommend his book "Differential Invariants."  Lie's style is a bit arcane and some modern rigor and nomenclature is missing, but the struggle is well worth it.  Many of his ideas have not been touched since his death.  
 A: The answer you present has no free parameters, while the ODE general solution should have $N$ free parameters.  But let's say the problem is to find any one non-zero special solution to the ODE. 
The easy way to get your answer is to start from the ansatz that the equation is solved by $y=kx^r$ for some $(k, r)$. Then 
$$
\frac{d^Ny}{dx^N} = k r^{\underline{N}}x^{r-N}
$$
where the falling power notation $r^{\underline{N}}$ means $r(r-1)(r-2)\cdots(r-N+1)$ with $N$ terms, so that $r^{\underline{1}} = r$ and $r^{\underline{r}}=r!$ if $r$ is a positive  integer.
We are solving
$$
y^{-M} \frac{d^Ny}{dx^N} = x^L \\
(kx^r)^{-M} k  r^{\underline{N}} x^{r-N}  = x^L \\
k^{1-M}r^{\underline{N}}x^{r-N-rM} = x^L
$$
Starting by equating the power of $x$ we have
$$
r-rM-N = L \\
r= \frac{L+N}{1-M} \\
k^{1-M} \left(\frac{L+N}{1-M}\right)^{\underline{N}} =1 \\
k = \left[\left(\frac{L+N}{1-M}\right)^{\underline{N}}\right]^\frac{1}{M-1} \\
y =  \left[\left(\frac{L+N}{1-M}\right)^{\underline{N}}\right]^\frac{1}{M-1}
x^{\frac{L+N}{1-M}}
$$
which matches your solution.
As to the general family of solutions, if your Lie theory method can solve that problem it is well worthwhile, because the equation is non-linear so the usual trick of taking the specific solution and adding a generic solution to a simpler equation does not work.
