Original Problem:
$$S=\displaystyle{\sum_{k=0}^{\infty}\binom{2k}{k}\frac{(-1)^k}{2^{2k}\left(4k+1\right)^{3}}}=\frac{1}{2}\int_{0}^{1}\frac{\ln^{2}x}{\sqrt{1+x^{4}}}dx$$
Variation I found:
$$\displaystyle{\sum_{k=0}^{\infty}\binom{2k}{k}\frac{1}{2^{2k}\left(4k+1\right)^{3}}}=\frac{1}{2}\int_{0}^{1}\frac{\ln^{2}x}{\sqrt{1-x^{4}}}dx$$
$$=\Gamma\left(\frac{1}{4}\right)^2\frac{16C+\pi^2}{128\sqrt{2\pi}}$$
where, $C$ is Catalan's Constant.
I am not able to Calculate the Integral:
$$I=\int_{0}^{1}\frac{\ln^{2}x}{\sqrt{1+x^{4}}}dx$$
Wolfram outputs a HyperGeometric Series as the Answer.
Best Answer
$$I=\int_0^1\frac{\ln^{2}x}{\sqrt{1+x^4}}dx=\frac12\int_0^\infty\frac{\ln^{2}x}{\sqrt{1+x^4}}dx=\frac1{128}\int_0^\infty\frac{x^{-\frac34}\ln^2x}{\sqrt{1+x}}dx$$ Denoting $$I_s=\frac1{128}\int_0^\infty\frac{x^{s-\frac34}}{\sqrt{1+x}}dx\,\,\Rightarrow\,\,I=\frac{\partial^2}{\partial s^2}\,\bigg|_{s=0}I_s$$ Then $$I_s\overset{t=\frac1{1+x}}{=}\frac1{128}\int_0^1t^{-s-\frac34}(1-t)^{s-\frac34}dt=\frac{\Gamma\big(\frac14-s\big)\Gamma\big(\frac14+s\big)}{128\sqrt\pi}$$ $$=\frac1{128\sqrt\pi}e^{\ln\Gamma\big(\frac14-s\big)+\ln\Gamma\big(\frac14+s\big)}$$ $$I=\frac{\partial^2}{\partial s^2}\,\bigg|_{s=0}I_s$$ $$=\frac{e^{\ln\Gamma\big(\frac14-s\big)+\ln\Gamma\big(\frac14+s\big)}}{128\sqrt\pi}\left(\Big(\psi\big(\frac14+s\big)-\psi\big(\frac14-s\big)\Big)^2+\psi^{(1)}\big(\frac14+s\big)+\psi^{(1)}\big(\frac14-s\big)\right)\,\bigg|_{s=0}$$ $$=\frac{\Gamma^2\big(\frac14\big)}{64\sqrt\pi}\psi^{(1)}\big(\frac14\big)=\frac{\Gamma^2\big(\frac14\big)(\pi^2+8G)}{64\sqrt\pi}=1.99282...$$
https://en.wikipedia.org/wiki/Trigamma_function