Another approach is obtained by introducing the Bessel functions to express
\begin{align}
f(x)&=\ln\left(\frac{\tan x-x}{x^3}\right)\\
&=\ln\left( \frac{\sin x-x\cos x}{x^3\cos x} \right)\\
&=\ln\left( \frac{J_{3/2}(x)}{x^2J_{-1/2}(x)} \right)
\end{align}
In a paper by Dickinson a uniformly convergent series expansion is given for the logarithm of a Bessel function:
\begin{equation}
\ln J_{\nu}(x)=\ln \left(\frac{(x / 2)^{\nu}}{\Gamma(\nu+1)}\right)+\sum_{n=1}^{\infty} \ln \left(1-\frac{x^2}{j_{\nu, n}}\right)
\end{equation}
when $\nu>-1$ and $0<x<j_{\nu ,1}$ ($j_{\nu,n}$ is the $n$-th zero of the Bessel function of order $\nu$). Then,
\begin{align}
\ln J_{\nu}(x)&=\ln \left(\frac{(x / 2)^{\nu}}{\Gamma(\nu+1)}\right)-\sum_{n=1}^{\infty}\sum_{k=0}^{\infty} \frac{1}{k+1}\frac{x^{2k+2}}{j_{\nu,n}^{2k+2}}\\
&=\ln \left(\frac{(x / 2)^{\nu}}{\Gamma(\nu+1)}\right)-\sum_{k=0}^{\infty} \frac{x^{2k+2}}{k+1}\sum_{n=1}^{\infty}\frac{1}{j_{\nu,n}^{2k+2}}
\end{align}
(all the summands in the convergent double series are positive and we may interchange the order of summation).
Thus
\begin{align}
f(x)&=-2\ln x+\ln J_{3/2}(x)-\ln J_{-1/2}(x)\\
&=-\ln3 +\sum_{k=0}^{\infty} \frac{x^{2k+2}}{k+1}\left(\sum_{n=1}^{\infty} \frac{1}{j_{-1/2,n}^{2k+2}}-\sum_{n=1}^{\infty}\frac{1}{j_{3/2,n}^{2k+2}} \right)
\end{align}
The sum of inverse even powers of zeros of Bessel functions appears in problems involving the diffusion equation and a broad literature exists on the topic. The quantity
\begin{equation}
\sigma(p,\nu)=\sum_{n=1}^\infty\frac{1}{j_{\nu,n}^{2p}}
\end{equation}
is discussed here. A recent analysis can be found in a paper by Jorge L. deLyra where a general formula for certain linear combinations of these sums is given ; it can be used to derive the formulas for $\sigma(p, \nu)$ by purely linear-algebraic means.
The first values for $\nu=3/2$ are
\begin{align}
&\sigma(1,3/2)=\frac{1}{10}\;;\;\sigma(2,3/2)=\frac{1}{350}\;;\;\sigma(3,3/2)=\frac{1}{7875}\\
&\sigma(4,3/2)=\frac{37}{6063750}\;;\;\sigma(5,3/2)=\frac{59}{197071875}
\end{align}
while $j_{-1/2,n}=(2n-1)\pi/2$. Then
\begin{align}
\sigma(k+1,-1/2)&=\sum_{n=1}^\infty\frac{1}{j_{-1/2,n}^{2k+2}}\\
&=\sum_{n=1}^\infty\left(\frac{2}{(2n-1)\pi} \right)^{2k+2}\\
&=\left( \frac{2}{\pi} \right)^{2k+2}(2^{2k+2}-1)\zeta(2k+2)\\
&=(2^{2k+2}-1)\frac{2^{2k+1}}{(2k+2)!}|B_{2k+2}|
\end{align}
Latter identity is obtained expressing the Zeta function at even values in terms of the Bernoulli numbers.
Finally,
\begin{align}
f(x)&=-\ln3 +\sum_{k=0}^{\infty} \frac{x^{2k+2}}{k+1}\left(\sigma(k+1,-1/2)-\sigma(k+1,3/2)\right)\\
&=-\ln3 +\sum_{k=0}^{\infty} \frac{x^{2k+2}}{k+1}\left((2^{2k+2}-1)\frac{2^{2k+1}}{(2k+2)!}|B_{2k+2}|-\sigma(k+1,3/2)\right)
\end{align}
which leads to the same first terms as in the expansion given by other users.
Series[Log[Tan[x] - x], {x, 0, 10}]in Mathematica.) Do you have any reason to believe that there is a better answer than this? $\endgroup$