I'm sure you've come across arXiv/1105.4312 --- this straightforward approach, numerical differentiation, is not necessarily less efficient than the route via the zeros. Alternatively, for that route, I think you'll find much of what you need in the work of Mark Coffey, in particular:
Relations and positivity results for the derivatives of the Riemann $\xi$ function (2004)
The $j$-th logarithmic derivative of the function $\xi(s)=(s/2)(s-1)\pi^{-s/2}\Gamma(s/2)\zeta(s)$ is given by
$$ \left[\frac{\xi'(s)}{\xi(s)}\right]^{(j)}=-\sum_{k=j+1}^{\infty}(k-1)(k-2)\cdots(k-j)\sigma_{k}s^{k-j-1} $$
where $\sigma_{j}=\sum_\rho \rho^{-j}$ is the sum of reciprocal powers of the complex zeros of the Riemann zeta function $\zeta(s)$. These coefficients are rigorously bounded by $O[(1+1/\pi)^j/j$, but actually decay much more rapidly. An efficient way to compute them numerically is given in the Appendix to the 2004 paper, in terms of a related coefficient $\eta_k$, related to $\sigma_k$ defined by $\sigma_k=(-1)^k\eta_{k-1}-(1-2^{-k})\zeta(k)+1$, $k\geq 2$.
