Here I present a new solution using the Burniston-Siewert method.
Let's start with the equation
\begin{equation}
\sin\left(x\right)-x\cos\left(x\right)=0
\end{equation}
which is equivalent, then we take the change of variable $x\rightarrow\textrm{arcsec}\left(z\right)+\pi k$ where $k\neq0$ we obtain the equation
\begin{equation}
\pi k+\textrm{arcsec}\left(z\right)-z\sqrt{1-1/z^{2}}=0
\end{equation}
Defining the function $f(z)=\pi k+\textrm{arcsec}\left(z\right)-z\sqrt{1-1/z^{2}}$ and applying the Burniston-Siewert method to solve transcendental
equations to this equation we obtain
\begin{equation}
x_{k}=\pi k+\textrm{arcsec}\left(\frac{\pi(2k+1)}{2}-m_{k}\right)
\end{equation}
where
\begin{equation}
m_{k}=\frac{1}{\pi}{\displaystyle \int\limits_{0}^{1}\theta_{k}\left(x\right)+\theta_{k+1}\left(x\right)\,dx\quad\textrm{and}\quad\theta_{k}\left(x\right)=\arctan\left(\frac{\textrm{arcsech}\left(x\right)-\sqrt{1-x^{2}}}{\pi k}\right)}
\end{equation}
or
\begin{equation}
m_{k}=\frac{1}{\pi}{\displaystyle \int\limits_{-1}^{1}\arctan\left(\frac{\textrm{arcsech}\left|x\right|-\sqrt{1-x^{2}}}{\pi/2\left(2k+1-\textrm{sgn}\left(x\right)\right)}\right)\,dx}
\end{equation}
Putting everything into a single expression
\begin{equation}
x_{k}=\pi k+\textrm{arcsec}\left(\frac{\pi(2k+1)}{2}-\frac{1}{\pi}{\displaystyle \int\limits_{-1}^{1}\arctan\left(\frac{\textrm{arcsech}\left|x\right|-\sqrt{1-1/x^{2}}}{\pi/2\left(2k+1-\textrm{sgn}\left(x\right)\right)}\right)\,dx}\right)
\end{equation}
These solutions are ordered with $k>0$, $x_{k}>0$ and $x_{k}<x_{k+1}$ (Negative solutions are obtained by multiplying by minus). From these expressions it can be seen that if $k\rightarrow\infty$ then $m_{k}\rightarrow0$ and then the roots approach to
\begin{equation}
\pi k+\textrm{arcsec}\left(\frac{\pi(2k+1)}{2}\right)
\end{equation}
Using the Laurent series of $\textrm{arcsec}\left(x\right)=\pi/2-1/x+\textrm{O}\left(x^{-2}\right)$ we obtain the most common approximation
\begin{equation}
x_{k}=\frac{\pi(2k+1)}{2}-\frac{2}{\pi(2k+1)}+\textrm{O}\left(k^{-2}\right)
\end{equation}
But, thanks to the representation obtained before we can find another better approximation. Using the fact that $\arctan\left(x\right)\sim x$ when $x$ as x approaches zero, we will have that
\begin{equation}
m_{k}=\frac{1}{\pi}{\displaystyle \int\limits_{-1}^{1}\frac{\textrm{arcsech}\left|x\right|-\sqrt{1-x^{2}}}{\pi/2\left(2k+1-\textrm{sgn}\left(x\right)\right)}\,dx}+\textrm{O}\left(k^{-3}\right)
\end{equation}
Integrating we obtain that
\begin{equation}
m_{k}=\frac{1}{4\pi}\cdot\frac{2k+1}{k\left(k+1\right)}+\textrm{O}\left(k^{-3}\right)
\end{equation}
Therefore
\begin{equation}
x_{k}=\pi k+\textrm{arcsec}\left(\left(2k+1\right)\left(\frac{4\pi^{2}k\left(k+1\right)-1}{8\pi k\left(k+1\right)}\right)\right)+\textrm{O}\left(k^{-5}\right)
\end{equation}
Finally I put a table comparing the real value $x_{k}$, the usual approximation $x_{k}^{1}$ and the improved approximation obtained before $x_{k}^{2}$
\begin{array}{c|ccc|cc}
\text{k} & x_{k} & x_{k}^{1} & x_{k}^{2} & \text{Rel. Err. } x_{k}^{1} & \text{Rel. Err. } x_{k}^{2} \\
\hline
1 & 4.493409457909064 & 4.500182389595496 & 4.492909579819157 & 0.2 & 0.01 \\
10 & 32.95638903982248 & 32.95640763543723 & 32.95638901560579 & 6.0 \times 10^{-5} & 7.35 \times 10^{-8} \\
100 & 315.7268944020431 & 315.7268944232251 & 315.7268944020428 & 6.71 \times 10^{-9} & 9.51 \times 10^{-14} \\
1000 & 3143.163131765756 & 3143.163131765777 & 3143.163131765756 & 6.83 \times 10^{-13} & 9.77 \times 10^{-20} \\
\end{array}
Ref:
-E. E. Burniston, C.E. Siewert. The use of Riemann problems in solving a class of transcendental equations. Mathematical Proceedings of the Cambridge Philosophical Society. 73. 111 - 118. (1973).
-Henrici P., Applied and computational complex analysis, Volume III. Wiley, New York. 183-191 (1986).
