Sum of the squares of the reciprocals of the fixed points of the tangent function
For the sake of at least having an answer, here's what's in J.M.'s link, in a nutshell anyway.
Note that $\tan(x)-x$ and $\sin(x)-x\cos(x)$ have the same zeros but the latter is holomorphic on the complex plane - the difference is that the latter has a triple root at $0$ whereas the former only has a double root, but our sum doesn't involve $x=0$ anyway. This means it affords a Weierstrass factorization. Some work can show that the coefficients of the expanded-out polynomials of the partial products in such a factorization will, in the limit, converge to those of the Taylor expansion of the function so long as we ignore the $e^{g(z)}$ and $E_n(z)$ factors. We can also show that the zeros are approximately $x_n\approx(n+\frac{1}{2})\pi$ asymptotically, and comparing this with the factorization for, say, $\sin$ tells us we don't need any $E_n$ factors. Now putting together the series expansions for sine and cosine gives
$$\left(x-\frac{1}{3!}x^3+\frac{1}{5!}x^5-\cdots\right)-x\left(1-\frac{1}{2!}x^2+\frac{1}{4!}x^4-\cdots\right)$$
$$=\frac{1}{3}x^3\left(1-\frac{1}{10}x^2+\cdots\right).$$
We ignore the $x^3$ and normalize so that $a_0=1,a_1=0,a_2=-1/10$. Now compute
$$\sum_{\tan(u)=u\ne0}u^{-2}=\left(\sum_{i}\frac{1}{\lambda_i}\right)^2-2\left(\sum_{i<j}\frac{1}{\lambda_i\lambda_j}\right)$$
$$=a_1^2-2a_2=\frac{1}{5}.$$
Note that $\tan(x)-x$ is an odd function and squaring takes out signs so by symmetry the above sum essentially double-counts every positive root. Divide by $2$ and our final answer is $1/10$.
Note that all sums-of-reciprocals of even powers of solutions to $\tan(x)=x$ can be evaluated in this way by using the Newton-Girard formulas.
The answer I got was not $\frac{1}{10}$ so it is likely that I have made a mistake somewhere.
$$x=\tan x$$
$$x^{2}\cos^{2}{x} =\sin^{2}{x}$$
$$x^{2}(1-\frac{x^{2}}{2!}+\frac{x^{4}}{4!}+\cdots)^{2}=(x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots)^{2}$$
$$x^{2}-x^{4}+\frac{x^{6}}{4}+\frac{x^{6}}{12}-\frac{x^{8}}{24}+\cdots=x^{2}-\frac{x^{4}}{3}+\frac{x^{6}}{36}+\frac{x^{6}}{60}+\cdots$$
$t=x^{2}$
$$t^{2}(\frac{2}{3}-\frac{13t}{45}+\cdots)=0$$
In a polynomial of the form $a_0+a_1x+\cdots+a_nx^{n}$ with roots $x_1,\ldots,x_n$, $\frac{a_0}{a_n}=x_{1}x_{2}\cdots$ and $\frac{a_2}{a_n}=-x_3x_4x_5-\cdots-x_1x_4-\cdots-x_1x_2x_5-\cdots-x_2x_4x_5-\cdots-\cdots$
If the root with $t=0$ is excluded, the sum of the reciprocals $=\frac{\frac{13}{45}}{\frac{2}{3}}=\frac{13}{30}$. $\tan{-x}=-\tan{x}$ therefore the reciprocal sums of positive and negative sums are equal, so for positive roots the sum is $\frac{13}{60}$.
I just saw a link here, so here's a solution of mine that dates back to 2008. Original version in this AoPS thread, which also includes the Weierstrass factorization solution. The solution here is modified from the original, most importantly by using a different contour so I don't have to worry about a principal value at a 3rd order pole.
Consider the function $$f(z)=\frac{\sin z}{z(\sin z-z\cos z)}=\frac{1}{z-z^2\cot z}$$ This function has simple poles at $z=x_n$ for nonzero roots $x_n=\tan x_n$ with residue $\frac{\sin x_n}{x_n}\cdot \frac1{x_n\sin x_n}=\frac1{x_n^2}$. Note that this residue is the same at the positive root $x_n$ and the negative root $-x_n$. It also has a third-order pole at $z=0$, where $f$ has Laurent series expansion \begin{align*}f(z) &= \frac1{z-z^2\frac{1-\frac12z^2+\frac1{24}z^4+\cdots}{z-\frac16z^3+\frac1{120}z^5+\cdots}}=\frac{z-\frac16z^3+\frac1{120}z^5+\cdots}{z^2-\frac16z^4+\frac1{120}z^6-z^2+\frac12z^4-\frac1{24}z^6+\cdots}\\ &= \frac{z-\frac16z^3+\frac1{120}z^5+\cdots}{\frac13z^4-\frac1{30}z^6+\cdots}=3z^{-3}\left(1-\frac16z^2+\cdots\right)\left(1+\frac1{10}z^2+\cdots\right)\\ f(z) &= 3z^{-3}-\frac15z^{-1}+\cdots\end{align*} for a residue of $-\frac15$.
Now, let $C(N,M)$ be the rectangular contour with corners at $N\pi+iM$, $-N\pi+iM$, $-N\pi-iM$, and $N\pi-iM$, for some large $M$ and large integer $N$. Define $I(N,M)=\int_{C(N,M)} f(z)\,dz$.
On the vertical segments, note that $\cot(\pm N\pi + iy)=\frac{\cos \pm N\pi\cosh y}{i\cos \pm N\pi\sinh y}=\frac1{i\tanh y}$. This is greater than $1$ in absolute value, so $|z^2\cot z-z|>|z|^2-|z|$ on those segments. Estimating $|z|\ge N\pi$, we get an integral of at most $\frac{2M}{N\pi(N\pi-1)}$ in absolute value on each vertical segment.
On the horizontal segments, $\cot(x+iy)$ tends to $\mp i$ uniformly as $y\to\pm\infty$. If $M$ is large enough that we're within $\epsilon$, this gives $|f(z)|\le \frac{1}{(1-\epsilon)^2|z|^2-|z|}$. Estimating $|z|\ge M$, this lead to an integral of at most $\frac{2N\pi}{M(M(1-\epsilon)^2-1)}$ in absolute value on each horizontal segment.
As long as neither $M$ nor $N$ gets too far ahead of the other - say, $M=\pi N$ - these both tend to zero as $M$ and $N$ go to $\infty$ together. We have shown that $\lim_{N\to\infty} I(N,\pi N)=0$.
Now, we calculate that integral with residues. Inside the contour $C(N,\pi N)$, there are $2N-1$ poles: the pole at zero, the first $N-1$ positive roots $x_1,x_2,\dots,x_{N-1}$ of $\tan x=x$, and their negatives $-x_1,-x_2,\dots,-x_{N-1}$. Therefore $$I(N,\pi N) = 2\pi i\left(-\frac15 + \sum_{k=1}^{N-1}\frac1{x_k^2}+\sum_{k=1}^{N-1}\frac1{(-x_k)^2}\right)$$ $$0 = 2\pi i\left(-\frac15 + 2\sum_{k=1}^{\infty}\frac1{x_k^2}\right)$$ $$\sum_{k=1}^{\infty}\frac1{x_k^2} = \frac1{10}$$
Connected : https://www.math.ucdavis.edu/~saito/publications/saito_rayleighfunc.pdf