If $I_n=\int_0^{\pi/4}\tan^n x dx$, then evaluate $\lim_{n\to\infty}\sum_{k=0}^nI_kI_{k+1}+I_kI_{k+3}+I_{k+1}I_{k+2}I_{k+3}$
$\color{brown}{\textbf{Theoretic data.}}$
Is known that
$$I_n = \int\limits_0^{\large\frac\pi4} \tan^n x\,\mathrm dx = \dfrac{\psi_0\left(\dfrac{n+3}4\right)-\psi_0\left(\dfrac{n+1}4\right)}4,\tag1$$
where $\psi_0(z)$ is digamma function,
$$\psi_0(z) = -\int\limits_0^\infty\dfrac{e^{-zt}}{1-e^{-t}}\,dt\tag2$$
(see also Wolfram Documentation).
So $$I_n = \frac14\int\limits_0^\infty\dfrac{1-e^{-\frac t2}}{1-e^{-t}}\,e^{\large-\frac {n+1}4t}\,dt = \frac18\int\limits_0^\infty\dfrac{e^{\large-\frac n4t}}{\cosh \dfrac t4}\,dt,$$ $$\color{green}{\mathbf{I_n = \frac12\int\limits_0^\infty\dfrac{e^{-nt}}{\cosh t}\,dt}}. \tag3$$
Taking in account OP, also is known $$I_0 = \frac\pi4,\quad I_1 = \frac12\ln2,\quad I_k+I_{k+2} = \dfrac1{k+1}.\tag4$$
$\color{brown}{\textbf{Summation of pairs.}}$
Applying immediate approach of OP, is possible to get
$$S_{20} = \sum\limits_{k=0}^{\infty}(I_kI_{k+1}+I_{k+1}I_{k+2})
= \sum\limits_{k=0}^{\infty}\dfrac{I_{k+1}}{k+1}
= \int\limits_0^{\large \frac\pi4}\sum\limits_{k=1}^{\infty}\dfrac{\tan^k
x}kdx
= -\int\limits_0^{\large \frac\pi4}\ln(1-\tan x)dx, $$
$$S_{20} = C - \frac\pi8\ln2\tag5$$
(see also Wolfram Alpha),
where C is Catalan's constant.
Then $$\sum\limits_{k=0}^{\infty}I_k(I_{k+1}+I_{k+3}) = \sum\limits_{k=0}^{\infty}\dfrac1{k+2}I_k = \sum\limits_{k=0}^{\infty}\dfrac1{k+2}\left(\dfrac1{k+1}-I_{k+2}\right)$$ $$= \sum\limits_{k=0}^{\infty}\left(\dfrac1{k+1}-\dfrac1{k+2}\right)-\sum\limits_{k=0}^{\infty}\dfrac1{k+2}I_{k+2} = 1+I_1-S_{20},$$ $$\color{green}{\boxed{\mathbf{S_2=\sum\limits_{k=0}^{\infty}(I_k I_{k+1}+I_k I_{k+3})= 1 + \frac{4+\pi}8\ln2-C}\approx 0.702\,806\,257.}}\tag6$$
$\color{brown}{\textbf{Integral presentation of product of external factors in the triples.}}$
Using $(3),$ one can write $$I_k I_{k+2} = \dfrac14\int\limits_0^\infty\int\limits_0^\infty \dfrac{e^{-k(x+y)}e^{-2y}}{\cosh x \cosh y}\,dy\,dx = \int\limits_0^\infty\int\limits_0^\infty \dfrac{e^{-(k+1)(x+y)}e^{4x}}{(1+e^{2x})(e^{2x+2y}+e^{2x})}\,dy\,dx.$$ Substitutions $$u=x+y,\quad t=e^{2x}$$ give $$I_k I_{k+2} = \dfrac12\int\limits_0^\infty e^{-(k+1)u} J_1(e^u)\,du,$$ where $$J_1(s) =\int\limits_1^{s^2}\dfrac{t}{(t+1)(t+s^2)}\,dt = \int\limits_1^{s^2}\dfrac1{s^2-1} \left(\dfrac {s^2}{t+s^2}-\dfrac{1}{t+1}\right)\,dt$$ $$ = \dfrac1{s^2-1}(s^2\ln(t+s^2)-\ln(t+1)\bigg|\phantom{\Big|}_1^{s^2}$$ $$= \dfrac{s^2\ln(2s^2) - (s^2+1)\ln(s^2+1) + \ln2}{s^2-1} = \dfrac{(s^2+1)\ln2+2s^2\ln s - (s^2+1)\ln(s^2+1)}{s^2-1}$$ $$= \dfrac{(s^2-1)\ln(s) - (s^2+1)\ln\dfrac{s^2+1}{2s}}{s^2-1} = \ln(s) - \dfrac{s^2+1}{s^2-1}\ln\dfrac{s^2+1}{2s}.$$
So
$$\color{green}{\mathbf{I_k I_{k+2}= \dfrac1{2(k+1)^2} - \frac12\int\limits_0^\infty\,e^{-(k+1)u}\coth u\ln\cosh u\,du.}} \tag7$$
Testing via formula (1) confirms formula(7)
$\color{brown}{\textbf{Summation of the triples productions}}$
Taking in account that $$\sum\limits_{k=2}^\infty \dfrac{t^k}{k^2} = \operatorname{Li}_2(t) - t$$ (see also Wolfram Alpha),
one can get $$\color{green}{\boxed{\mathbf{S_{30} = \sum\limits_{k=2}^\infty \dfrac{I_k}{k^2} =\int\limits_0^{\large\frac\pi4}(\operatorname{Li}_2(\tan t)-\tan t)\, dt} \approx 0.090\,057\,376}} \tag{8a}$$ (see also Wolfram Alpha numeric calculations).
The same result can be obtained via summation: $$S_{30} = \sum\limits_{k=2}^\infty \dfrac{I_k}{k^2} = \frac14\sum\limits_{k=2}^\infty \frac1{k^2}\left(\psi_0\left(\dfrac{k+3}4\right)-\psi_0\left(\dfrac{k+1}4\right)\right)$$ $$= \frac14\sum\limits_{k=3}^\infty \left(\frac1{(k-1)^2}-\frac1{(k+1)^2}\right)\psi_0\left(\dfrac{k+2}4\right) - \dfrac1{16}\psi\left(\frac34\right) - \dfrac1{36}\psi_0(1)$$ $$= \sum\limits_{k=3}^\infty \frac k{(k^2-1)^2}\psi_0\left(\dfrac{k+2}4\right) + \dfrac1{16}\left(\gamma+3\ln2-\dfrac\pi2\right) + \dfrac\gamma{36},$$ where $\gamma$ is Euler-Maclaurin constant. Taking in account that $$\dfrac{13\gamma}{144}+\dfrac{3\ln2}{16}-\dfrac\pi{32}\approx0.083\,900\,073\,4$$ (see also Wolfram Alpha) and $$\sum\limits_{k=3}^\infty \frac k{(k^2-1)^2}\psi_0\left(\dfrac{k+2}4\right) \approx0.006\,157\,302\,6$$ (see also Wolfram Alpha) one can get $$\color{green}{\boxed{\mathbf{S_{30}=\dfrac{13\gamma}{144}+\dfrac{3\ln2}{16}-\dfrac\pi{32} +\sum\limits_{k=3}^\infty \frac k{(k^2-1)^2}\psi_0\left(\dfrac{k+2}4\right)}\approx0.090\,057\,376.}}\tag{8b}$$
The second sum is $$S_{31} = \sum\limits_{k=1}^\infty I_{k+1} \int\limits_0^\infty\,e^{-(k+1)u}\coth u\ln\cosh u\,du = \sum\limits_{k=1}^\infty \int\limits_0^\infty\,\int\limits_0^\infty\, e^{-(k+1)(u+v)}\dfrac{\coth u}{\cosh v}\ln\cosh u\,dv\,du$$ $$ = \int\limits_0^\infty\,\int\limits_0^\infty\, \dfrac{e^{-2(u+v)}}{1-e^{-(u+v)}} \dfrac{\coth u}{\cosh v}\ln\cosh u\,dv\,du$$ $$ = \int\limits_0^\infty\,\int\limits_0^\infty\, \dfrac{e^{-2u-3v}}{1-e^{-(u+v)}} \cdot \dfrac{\coth u}{1+e^{-2v}}\ln\cosh u\,dv\,du = \int\limits_0^\infty J_2(e^{-u})\coth u\ln\cosh u\,du,$$ where $$J_2(x) = \int\limits_0^1 \dfrac{x^2t^2\,dt}{(1-xt)(t^2+1)} = \dfrac {x^2}{x^2+1}\int\limits_0^1 \left(\dfrac{1}{1-xt} - \dfrac{xt+1}{t^2+1}\right) \,dt$$ $$ = -\dfrac x{x^2+1}\left(\ln(1-xt) + \dfrac{x^2}2\ln(t^2+1) + x\arctan t\right) \bigg|\phantom{\Big|}_0^1 = \dfrac18 \dfrac{2x}{x^2+1} \left(-4\ln(1-x) - 2x^2\ln2 - \pi x\right).$$
Then $$\color{green}{\boxed{\mathbf{S_{31} = \int\limits_0^\infty \dfrac{4e^{2u}\ln(1-e^{-u}) - \pi e^u -2\ln2}{8e^{3u}\sinh u}\ln\cosh u\,du}\approx 0.025\,054\,224}}\tag9$$ (see also Wolfram Alpha numeric calculations)
$\color{brown}{\textbf{Final summations.}}$
$$\boxed{S_3 = \sum\limits_{k=1}^\infty I_k I_{k+1} I_{k+2} = \frac12(S_{30} - S_{31})\approx \frac12(0.090\,057\,376 - 0.025\,054\,224) = 0.032\,501\,576,}$$ $$\color{green}{\boxed{\phantom{\Bigg|}S=S_2+S_3\approx 0.702\,806\,257+0.032\,501\,576 = \color{brown}{\mathbf{0.735\,307\,833}.\quad}}}$$
\begin{aligned} I_{n} &=\int \tan ^{n} x d x \\ &=\int \tan ^{n-2} x \tan ^{2} x d x \\ &=\int \tan ^{n-2} x\left(\sec ^{2} x-1\right) d x \\ &=\int \tan ^{n-2} x \sec ^{2} x d x-I_{n-2} \end{aligned}
Use Integration by parts: \begin{aligned} I_{n} &=\tan ^{n-1} x-\int(n-2) \tan ^{n-2} x \sec ^{2} x d x \\ &=\tan ^{n-1} x-(n-2)\left(I_{n}+I_{n-2}\right)-I_{n-2} \end{aligned}
Which gives $$ I_{n}=\frac{1}{n-1} \tan ^{n-1} x-I_{n-2} $$
Iteratively,
$$ I_{2 k+1}=(-1)^{k} I_{1}+\sum_{p=1}^{k} \frac{(-1)^{p+k}}{2 p} \tan ^{2 p} x $$ and $$ I_{2 k}=(-1)^{k} I_{0}+\sum_{p=1}^{k} \frac{(-1)^{p+k}}{2 p-1} \tan ^{2 p-1} x $$
We need the boundary values of $I_0$ and $I_1$
$$ \begin{aligned} I_{0} &=\int \tan ^{0} x d x \\ &=x \end{aligned} $$ and $$ \begin{aligned} I_{1} &=\int \tan ^{1} x d x \\ &=\int \tan x d x \\ &=-\ln \cos x \\ &=\ln \sec x \end{aligned} $$
Therefore, the final expression for any $I_n$ would be $$ I_{2 k+1}=(-1)^{k} \ln \sec x+\sum_{p=1}^{k} \frac{(-1)^{p+k}}{2 p} \tan ^{2 p} x $$ and
$$ I_{2 k}=(-1)^{k} x+\sum_{p=1}^{k} \frac{(-1)^{p+k}}{2 p-1} \tan ^{2 p-1} x $$
Now, $\tan^n(\frac{pi}{4}) = 1$, we end up with partial fractions.
$$ I_{2 k+1}=(-1)^{k} \ln \sqrt{2}+\sum_{p=1}^{k} \frac{(-1)^{p+k}}{2 p} $$ and
$$ I_{2 k}=(-1)^{k} \frac{\pi}{4}+\sum_{p=1}^{k} \frac{(-1)^{p+k}}{2 p-1} $$
I am unable to simplify this further. Please feel free to edit it.