Is there an "elegant" non-recursive formula for these coefficients? Also, how can one get proofs of these patterns?
Let $\beta_n$ denote the flag $h$-vector (as defined in EC1, Section 3.13) of the partition lattice $\Pi_n$ (EC1, Example 3.10.4). Then $$ \mathrm{mag}_{n,{n-1\choose 2}-j} = \sum_S \beta_n(S), $$ where $S$ ranges over all subsets of $\lbrace 1,2,\dots,n-2\rbrace$ whose elements sum to $j$. An explicit formula for $\beta_n(S)$ is given by $$ \beta_n(S) = \sum_{T\subseteq S} (-1)^{|S-T|} \alpha_n(T), $$ where if the elements of $T$ are $t_1<\cdots < t_k$, then $$ \alpha_n(T) = S(n,n-t_1)S(n-t_1,n-t_2) S(n-t_2,n-t_3)\cdots S(n-t_{k-1},n-t_k). $$ Here $S(m,j)$ denotes a Stirling number of the second kind.
Addendum. A combinatorial description of the mag numbers is somewhat complicated. Consider all ways to start with the $n$ sets $\lbrace 1 \rbrace,\dots, \lbrace n \rbrace$. At each step we take two of our sets and replace them by their union. After $n-1$ steps we will have the single set $\lbrace 1,2,\dots,n \rbrace$. An example for $n=6$ is (writing a set like $\lbrace 2,3,5\rbrace$ as 235) 1-2-3-4-5-6, 1-2-36-4-5, 14-36-2-5, 14-356-2, 14356-2, 123456. At the $i$th step suppose we take the union of two sets $S$ and $T$. Let $a_i$ be the least integer $j$ such that $j$ belongs to one of the sets $S$ or $T$, and some number smaller than $j$ belongs to the other set. For the example above we get $(a_1,\dots,a_5)=(6,4,5,3,2)$. If $\nu$ denotes this merging process, then let $f(\nu) = \sum i$, summed over all $i$ for which $a_i>a_{i+1}$. For the above example, $f(\nu) = 1+3+4=8$. (The number $f(\nu)$ is called the major index of the sequence $(a_1,\dots,a_{n-1})$.) Then $\mathrm{mag}_{n,{n-1\choose 2}-j}$ is the number of merging processes $\nu$ for which $f(\nu)=j$. This might look completely contrived to the uninitiated, but it is very natural within the theory of flag $h$-vectors.
Mike, possibly you know this all, but for my own pleasure (and for the reader who's not yet much familiar with this all) I can give an explicite description for the mag-numbers. However, this is simply an explication of the terms of some series/arrays which are in fact recursive, but the recursion is so flat that we can resolve it without too much hazzle to direct references on factorials/binomials and the log of the fixpoint only. I employ the notation of (operator-)matrices which are known as Bell-/Carleman-matrices.
The text became too long for the answer-field here, so I link to a pdf-file on my homepage.
(If you don't like pdf there is also a html-version, but automatically generated by word and not perfect formatted)
Since I describe this using the known procedure of diagonalization of a triangular matrix some of your questions concerning the structure of coefficients may be answered or possibly a rigorous answer lays on the hand.
(P.s.: if it is more convenient for the MO-readers I could upload or possibly reformat the text for mathjax, but the latter would be much unwanted work...)
[update]: updated the pdf-file for readability
Bell polynomials
The first step to solving the problem is to deal with the recursion combinatorially. Note the Bell polynomial in $$a_n = \frac{B_n(1! a_1, 2! a_2, ..., (n-1)! a_{n-1}, 0)}{n!(L^{n-1} - 1)}.$$
Bell polynomials are generalizations of the combinatorial structure set partitions or Bell numbers. The Faà di Bruno's formula in terms of the Bell polynomials provides a way to express the Taylor's series of a composite function $f(g(x))$.
From Wikipedia,
$${d^n \over dx^n} f(g(x)) = \sum_{k=1}^n f^{(k)}(g(x))\cdot B_{n,k}\left(g'(x),g''(x),\dots,g^{(n-k+1)}(x)\right).$$
The formula has a "combinatorial" form:
$${d^n \over dx^n} f(g(x))=(f\circ g)^{(n)}(x)=\sum_{\pi\in\Pi} f^{(\left|\pi\right|)}(g(x))\cdot\prod_{B\in\pi}g^{(\left|B\right|)}(x)$$
where
π runs through the set Π of all partitions of the set { 1, ..., ''n'' },
$B ∈ π$ means the variable $B$ runs through the list of all of the "blocks" of the partition π, and
$|A|$ denotes the cardinality of the set $A$ (so that $|π|$ is the number of blocks in the partition $π$ and $|B|$ is the size of the block $B$).
The appearance of Bell polynomials in this problem is not surprising as the problem is an example of fractional iteration where $f(g(x))=f(f^{m-1}(x)).$
Combinatorial Structures
Moving the recursion from the analytic to the combinatoric, given that the set partitions are the combinatorial form of the derivatives of composite functions, what is the combinatorial form associated with the derivatives of $D(f(g(x)))=D^nf^m(x)$, iterated functions? Because there are now an arbitrary number of levels of recursion, there are also an arbitrary number of recursive set partitionings. This is equivalent to the combinatorial problem of the total partitions of $n$ or Schroeder's Fourth Problem.
Analytic
Without loss of generalization assume that $f(0)=g(0)=0$. The following proof removes with restriction that $a_1=1$. It will be seen that this approach handles both the cases of Schroeder's equation and Abel's equation.
The First Derivative
The first derivative of a function at its fixed point $Df(0)=f_1$ is often represented by $\lambda$ and referred to as the multiplier or the Lyapunov characteristic number; its logarithm is known as the Lyapunov exponent. Let $g(z)=f^{m-1}(z)$, then
$$Df(g(z))=f'(g(z))g'(z)$$
$$=f'(f^{m-1}(z))Df^{m-1}(z)$$
$$=\prod^{m-1}_{k_1=0}f'(f^{m-k_1-1}(z))$$
$$Df^m(0)=f'(0)^m$$
$$=f_1^m = \lambda^m$$
The Second Derivative
$$D^2f(g(z))=f''(g(z))g'(z)^2+f'(g(z))g''(z)$$ $$=f''(f^{m-1}(z))(Df^{m-1}(z))^2+f'(f^{m-1}(z))D^2f^{m-1}(z)$$
Setting $g(z) = f^{m-1}(z)$ results in $$D^2f^m(0)= f_2 \lambda^{2m-2}+\lambda D^2f^{m-1}(0)$$ When $\lambda \neq 0$, a recurrence equation is formed that is solved as a summation. Note that $|\lambda|\neq1$ is consistent with Schroder's equation, while $\lambda=1$ is consistent with Abel's equation. $$ D^2f^m(0)=f_2\lambda^{2m-2}+\lambda D^2f^{m-1}(0)$$ $$ =\lambda^0f_2 \lambda^{2m-2} +\lambda^1f_2 \lambda^{2m-4} +\cdots +\lambda^{m-2}f_2 \lambda^2 +\lambda^{m-1}f_2 \lambda^0 $$ $$ =f_2\sum_{k_1=0}^{m-1}\lambda^{2m-k_1-2} $$
The Third Derivative
Continuing on with the third derivative,$$ D^3f(g(z))=f'''(g(z))g'(z)^3+3f''(g(z))g'(z)g''(z)+f'(g(z))g'''(z) =f'''(f^{m-1}(z))(Df^{m-1}(z))^3 +3f''(f^{m-1}(z))Df^{m-1}(z)D^2f^{m-1}(z) +f'(f^{m-1}(z))D^3f^{m-1}(z) $$
$$ D^3f^m(0)=f_3\lambda^{3m-3}+3 f_2^2\sum_{k_1=0}^{m-1}\lambda^{3m-k_1-5}+\lambda D^3f^{m-1}(0)$$ $$ =f_3\sum_{k_1=0}^{m-1}\lambda^{3m-2k_1-3} +3f_2^2 \sum_{k_1=0}^{m-1} \sum_{k_2=0}^{m-k_1-2} \lambda^{3m-2k_1-k_2-5} $$ Note that the index $k_1$ from the second derivative is renamed $k_2$ in the final summation of the third derivative. A certain amount of renumbering is unavoidable in order to use a simple index scheme.
The Higher Derivatives
The fourth derivative is comprised of five summations.
$D^n f^m(0) = \sum \frac{n!(D^k f)(0)}{k_1! \cdots k_{n-1}!} $ $\left(\frac{Df^{m-1}(0)}{1!}\right)^{k_1} \cdots $ $\left(\frac{D^n f^{m-1}(0)}{(n-1)!}\right)^{k_{n-1}} $ $ + \lambda D^n f^{m-1}(0)$
Computation
SchroederSummations.nb and Iterate.m are Mathematica notebooks which are an attempt to systematically generate the combinatorial structure of unlabeled total partitions and then compute the analytic form of each combinatorial structure. Tallying the integer coefficients of the derivatives gives the sequence 1,1,4,26,236,2752,39208,660032,$\ldots$ which is the structure total partitions. The number of summations is given by the structure unlabeled total partitions; 1, 1, 2, 5, 12, 33, 90, 261,$\ldots$.
So the eighth derivative with 660032 terms, 261 with symmetry, can be directly computed without the lesser derivatives. Faà di Bruno's formula is indexed by the set partitions, but the derivatives of an iterated function are indexed by total partitions. Each total partition maps into a summation.
Note: I did spend some time trying to find a combinatorial structure associated with $a_n$. I tried setting L to a few reasonable values; -1, -1/2, 1/2 and 1 to see if I could generates a combinatorial sequence in the OEIS, but I didn't find any matching sequences.