Derivation of asymptotic solution of $\tan(x) = x$. An equation that seems to come up everywhere is the transcendental $\tan(x) = x$. Normally when it comes up you content yourself with a numerical solution usually using Newton's method. However, browsing today I found an asymptotic formula for the positive roots $x$:
$x = q - q^{-1} - \frac23 q^{-3} + \cdots$
with $q = (n + 1/2) \pi$ for positive integers $n$. For instance here: http://mathworld.wolfram.com/TancFunction.html, and here: http://mathforum.org/kb/message.jspa?messageID=7014308 found from a comment here: Solution of tanx = x?.
The Mathworld article says that you can derive this formula using series reversion, however I'm having difficulty figuring out exactly how to do it.
Any help with a derivation would be much appreciated.
 A: You may be interested in N. G. de Bruijn's book Asymptotic Methods in Analysis, which treats the equation $\cot x = x$.  What follows is essentially a minor modification of that section in the book.
The central tool we will use is the Lagrange inversion formula.  The formula given in de Bruijn differs slightly from the one given on the wiki page so I'll reproduce it here.

Lagrange Inversion Formula.
Let the function $f(z)$ be analytic in some neighborhood of the point $z=0$ of the complex plane.  Assuming that $f(0) \neq 0$, we consider the equation $$w = z/f(z),$$ where $z$ is the unknown.  Then there exist positive numbers $a$ and $b$ such that for $|w| < a$ the equation has just one solution in the domain $|z| < b$, and this solution is an analytic function of $w$: $$z = \sum_{k=1}^{\infty} c_k w^k \hspace{1cm} (|w| < a),$$ where the coefficients $c_k$ are given by $$c_k = \frac{1}{k!} \left\{\left(\frac{d}{dz}\right)^{k-1} (f(z))^k\right\}_{z=0}.$$

Essentially what this says is that we can solve the equation $w = z/f(z)$ for $z$ as a power series in $w$ when $|w|$ and $|z|$ are small enough.
Okay, on to the problem.  We wish to solve the equation $$\tan x = x.$$  As with many asymptotics problems, we need a foothold to get ourselves going.  Take a look at the graphs of $\tan x$ and $x$:

We see that in each interval $\left(\pi n - \frac{\pi}{2}, \pi n + \frac{\pi}{2}\right)$ there is exactly one solution $x_n$ (i.e. $\tan x_n = x_n$), and, when $n$ is large, $x_n$ is approximately $\pi n + \frac{\pi}{2}$.  But how do we show this second part?
Since $\tan$ is $\pi$-periodic we have
$$\tan\left(\pi n + \frac{\pi}{2} - x_n\right) = \tan\left(\frac{\pi}{2} - x_n\right)$$
$$\hspace{2.4 cm} = \frac{1}{\tan x_n}$$
$$\hspace{2.6 cm} = \frac{1}{x_n} \to 0$$
as $n \to \infty$, where the second-to-last equality follows from the identites $$\sin\left(\frac{\pi}{2} - \theta\right) = \cos \theta,$$ $$\cos\left(\frac{\pi}{2} - \theta\right) = \sin \theta.$$
Since $-\frac{\pi}{2} < \pi n + \frac{\pi}{2} - x_n < \frac{\pi}{2}$ and since $\tan$ is continuous in this interval we have $\pi n + \frac{\pi}{2} - x_n \to 0$ as $n \to \infty$.  Thus we have shown that $x_n$ is approximately $\pi n + \frac{\pi}{2}$ for large $n$.
Now we begin the process of putting the equation $\tan x = x$ into the form required by the Lagrange inversion formula.  Set $$z = \pi n + \frac{\pi}{2} - x$$ and $$w = \left(\pi n + \frac{\pi}{2}\right)^{-1}.$$  Note that we do this because when $|w|$ is small (i.e. when $n$ is large) we may take $|z|$ small enough such that there will be only one $x$ (in the sense that $x = \pi n + \frac{\pi}{2} - z$) which satisfies $\tan x = x$.  Plugging $x = w^{-1} - z$ into the equation $\tan x = x$ yields, after some simplifications along the lines of those already discussed, $$\cot z = w^{-1} - z,$$ which rearranges to $$w = \frac{\sin z}{\cos z + z\sin z} = z/f(z),$$ where $$f(z) = \frac{z(\cos z + z\sin z)}{\sin z}.$$  Here note that $f(0) = 1$ and that $f$ is analytic at $z = 0$.  We have just satisfied the requirements of the inversion formula, so we may conclude that we can solve $w = z/f(z)$ for $z$ as a power series in $w$ in the form given earlier in the post.
We have $c_1 = 1$ and, since $f$ is even, it can be shown that $c_{2k} = 0$ for all $k$.  Calculating the first few coefficients in Mathematica gives $$z = w + \frac{2}{3}w^3 + \frac{13}{15}w^5 + \frac{146}{105}w^7 + \frac{781}{315}w^9 + \frac{16328}{3465}w^{11} + \cdots.$$  Substituting this into $x = w^{-1} - z$ and using $w = \left(\pi n + \frac{\pi}{2}\right)^{-1}$ gives the desired series for $x_n$ when $n$ is large enough: $$x_n = \pi n + \frac{\pi}{2} - \left(\pi n + \frac{\pi}{2}\right)^{-1} - \frac{2}{3}\left(\pi n + \frac{\pi}{2}\right)^{-3} - \frac{13}{15}\left(\pi n + \frac{\pi}{2}\right)^{-5} - \frac{146}{105}\left(\pi n + \frac{\pi}{2}\right)^{-7} - \frac{781}{315}\left(\pi n + \frac{\pi}{2}\right)^{-9} - \frac{16328}{3465}\left(\pi n + \frac{\pi}{2}\right)^{-11} + \cdots$$
