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}$$