Writing $\tan^{-1}\cot(\pi x)$ is in fact just a very fancy/silly way to write the sawtooth function.
Let $B_1(x)=s(x)=\{x\}-\frac{1}{2}$ denote the sawtooth function. Then note that $$\tan^{-1}\cot(\pi x)=\pi^2 s(x)^2.$$ Their identity is then $$\zeta(s)=\frac{s}{s-1}-\frac{1}{2}+\frac{s}{8}-\frac{s(s+1)}{2}\int_{1}^{\infty} \left(\{u\}-\frac{1}{2}\right)^2 u^{-s-2}du,$$ however, as you have noted, this isn't correct either. So what is the real identity lurking here, and how do we prove it? In fact, we have that $$\zeta(s)=\frac{s}{s-1}-\frac{1}{2}-\frac{s(s+1)}{2}\int_{1}^{\infty}\left(\{u\}^{2}-\{u\}\right)u^{-s-2}du $$ and this follows from integration by parts twice.
Writing $$\zeta(s)=\sum_{n=1}^{\infty}n^{-s}=\int_{1}^{\infty}x^{-s}d\lfloor x\rfloor,$$ and applying integration by parts, we have $$\zeta(s)=\frac{s}{s-1}-s\int_{1}^{\infty}\{u\}u^{-s-1}du.$$ Now, we could proceed from here, but instead we'll rewrite this in terms of the first periodic bernoulli polynomial $s(x)$ to obtain $$\zeta(s)=\frac{s+1}{2\left(s-1\right)}-s\int_{1}^{\infty}B_{1}(\{u\})u^{-s-1}du.$$
Applying integration by parts again, $$\int_{1}^{\infty}B_{1}(\{u\})u^{-s-1}du=u^{-s-1}\int_{1}^{u}B_{1}(\{u\})dt\biggr|_{1}^{\infty}+(s-1)\int_{1}^{\infty}u^{-s-2}\left(\int_{1}^{u}B_{1}(\{u\})dt\right)du.$$ Now, as $$\int_{0}^{x}B_{1}(\{u\})du=\frac{B_{2}\left(\{x\}\right)-B_{2}(0)}{2}=\frac{1}{2}\left\{ x\right\} ^{2}-\frac{1}{2}\left\{ x\right\},$$ we arrive at the desired identity. We could continue in this way, and rewrite the identity as $$\zeta(s)=\frac{s}{s-1}-\frac{1}{2}+\frac{s}{12}-\frac{s(s+1)}{2}\int_{1}^{\infty}B_{2}(\{u\})u^{-s-2}du,$$ applying integration by parts again, using the fact that $$\int_0^x B_n(\{t\})dt = \frac{B_{n+1}({x})-B_{n+1}(0)}{n+1},$$ and this yields a method of extending $\zeta(s)$ to the half plane $\text{Re}(s)>-n$ where the formula involves the first $n+1$ Bernoulli numbers. In fact, some combinatorial analysis can turn this into an alternate proof that $$\zeta(-n)=-\frac{B_{n+1}}{n+1},$$ and hence if combined with the functional equation for the zeta function we obtain an alternate proof for the value of $\zeta(2k)$.
Your first sum equals
$$
\sum_{i,j}\frac1{i^nj^n}-\sum_{i=j}\frac1{i^nj^n}=\zeta(n)^2-\zeta(2n).$$
The same trick works for your the second sum.
Best Answer
In this MSE question the OP derives to a very similar expression $$x\cot x = 1-2\sum_{n=1}^{\infty}\zeta(2n)\frac{x^{2n}}{\pi^{2n}}.$$ To get to this result he starts with the Weiertrass product of the sine. So, how can we get from this expression to yours? First of all we can bring the $1$ as well as the $-2$ to the left hand side. $$\frac{1-x\cot x}{2} = \sum_{n=1}^{\infty}\zeta(2n)\frac{x^{2n}}{\pi^{2n}}.$$ Now substitute $x\to\pi x$. This yields $$\frac{1-\pi x\cot \pi x}{2} =\sum_{n=1}^{\infty}\zeta(2n)\frac{(\pi x)^{2n}}{\pi^{2n}}=\sum_{n=1}^{\infty}\zeta(2n)x^{2n}.$$
Start from $n=0$ instead from $n=1$, by writing every $n$ as $n+1$. $$\sum_{n=1}^{\infty}\zeta(2n)x^{2n}=\sum_{n=0}^{\infty}\zeta(2(n+1))x^{2(n+1)} =x^2\sum_{n=0}^{\infty}\zeta(2n+2)x^{2n}.$$ Remember what the sum was equal to and divide both sides by $x^2$ $${\frac{1-\pi x\cot \pi x}{2x^2} = \sum_{n=0}^{\infty}\zeta(2n+2)x^{2n}}.$$