Let $d \ge 1$ and $n \ge d+1$ be integers. We consider a following multiple sum: \begin{equation} {\mathfrak S}^{(d)}(n):= \sum\limits_{1 \le n_1 < n_2 < \cdots < n_d \le n-1} \prod\limits_{\xi=1}^{d+1} \frac{1}{(2 (n_\xi-n_{\xi-1})+1)!} \end{equation} subject to $n_0=0$ and $n_{d+1}=n$. By using the recursion relation \begin{equation} {\mathfrak S}^{(d+1)}(n) = \sum\limits_{n_1=1}^{n-1} \frac{1}{(2 n_1+1)!} \cdot {\mathfrak S}^{(d)}(n-n_1) \end{equation} and then by using the binomial expansion formula we have found the following results: \begin{eqnarray} &&{\mathfrak S}^{(d)}(n) = \\ \frac{1}{(2n+1)!} \quad && \quad \mbox{if $d=0$}\\ \frac{2^{2n+2}}{2 (2n+2)!} - \frac{2}{(2n+1)!} \quad &&\quad \mbox{if $d=1$}\\ \frac{3^{2n+3}-3}{4(2n+3)!}-\frac{3\cdot 2^{2n+2}}{2(2n+2)!}+\frac{3}{(2n+1)!}\quad && \quad \mbox{if $d=2$}\\ \frac{4^{2n+4}-4\cdot 2^{2n+4}}{8(2n+4)!}-\frac{4 \cdot 3^{2n+3}-12}{4 (2n+3)!}+\frac{6\cdot 2^{2n+2}}{2(2n+2)!}-\frac{4}{(2n+1)!} \quad && \quad \mbox{if $d=3$}\\ \frac{5^{2n+5}-5 \cdot 3^{2n+5}+10}{16(2n+5)!}-\frac{5\cdot 4^{2n+4}-20 \cdot 2^{2n+4}}{8(2n+4)!}+\frac{10 \cdot 3^{2n+3}-30}{4(2n+3)!}-\frac{10 \cdot 2^{2n+2}}{2(2n+2)!}+\frac{5}{(2n+1)!}\quad && \quad \mbox{if $d=4$}\\ \frac{6^{2n+6}-6 \cdot 4^{2n+6}+15 \cdot 2^{2n+6}}{32 (2n+6)!}-\frac{6 \cdot 5^{2n+5}-30 \cdot 3^{2n+5}+60}{16 (2n+5)!}+\frac{15 \cdot 4^{2n+4}-60 \cdot 2^{2n+4}}{8(2n+4)!}-\frac{20 \cdot 3^{2n+3}-60}{4(2n+3)!}+\frac{15 \cdot 2^{2n+2}}{2(2n+2)!}-\frac{6}{(2n+1)!} \quad && \quad \mbox{if $d=5$} \end{eqnarray} Now, by using the above results and the expansion: \begin{eqnarray} \frac{z}{\sinh(z)}&=& 1+ \sum\limits_{n=1}^\infty \left(\sum\limits_{d=0}^{n-1} (-1)^{d+1} \cdot {\mathfrak S}^{(d)}(n)\right) z^{2 n}\\ &=&1-\frac{z^2}{6}+\frac{7 z^4}{360}-\frac{31 z^6}{15120}+\frac{127 z^8}{604800}-\frac{73 z^{10}}{3421440}+\frac{1414477 z^{12}}{653837184000}+O(z^{14}) \end{eqnarray} we have found the expansion coefficients up to the twelve-th order.
Now the question is twofold. Firstly, how does the formula for our multiple sum look like for generic values of $d$ and secondly, how else can we derive a Taylor expansion of the function $z/sinh(z)$?
It is not difficult to see a pattern in the multiple sums in the formulation of the question. Therefore the generic result reads: \begin{eqnarray} {\mathfrak S}^{(d)}(n) = \sum\limits_{j=1}^{d+1} \sum\limits_{p=0}^{\lfloor j/2 \rfloor}\frac{(-1)^{d+1-j-p}}{2^{j-1}(2n+j)!} \binom{d+1}{j} (j-2 p)^{2n+j} \binom{j}{p} \end{eqnarray} for $d\ge 0$ and $n \ge d+1$. By summing the result over $d$ (with the appropriate sign) we get the Taylor expansion being sought for: \begin{eqnarray} \frac{z}{\sinh(z)} = 1+\sum\limits_{n=1}^\infty (-1)^n {\mathcal a}_n z^{2 n} \end{eqnarray} where \begin{equation} {\mathcal a}_n := \left( \sum\limits_{j=1}^n \sum\limits_{p=0}^{\lfloor j/2 \rfloor } \frac{(-1)^{j+p}}{2^{j-1} (2n+j)!} \binom{n+1}{j+1} (j-2 p)^{2n+j} \binom{j}{p}\right) \end{equation} Interestingly enough we have checked numerically that the limit of consecutive ratios converges as follows: \begin{equation} \lim\limits_{n\rightarrow \infty} \frac{a_{n+1}}{a_n} = \frac{1}{\pi^2} \end{equation} as it should be.