If you choose to evaluate this using a contour it will not just be similar, but identical to the contour for $\cot(\pi z)$ with a variable change.
It is literally identical, with only a variable change, and you can also choose to do this variable change before or after the integration itself. What you are asking for in the question is doing the variable change before the contour integration, what follows below is doing the variable change after:
By contour integration we have that $$\pi z\cot(\pi z)=1+2z^2\sum_{n=1}^\infty \frac{1}{z^2-n^2}.$$ Let $z=\frac{i}{2}$, then this is $$ \frac{\pi i}{2}\cot\left(\frac{\pi i}{2}\right)=1+2\sum_{n=1}^\infty \frac{1}{4n^2+1}.$$ Since $\coth(z)=i\cot(iz)$ we conclude that $$\sum_{n=0}^\infty \frac{1}{4n^2+1}=\frac{1}{2}+\frac{\pi}{4}\coth\left(\frac{\pi}{2}\right).$$ as desired.
If you want to do the variable change before the contour integration, just look at $z\rightarrow \frac{iz}{2}$. Then we are integrating the function $\pi z \coth \left( \frac{\pi z}{2}\right)$, which really the same function as $\pi z\cot\left(\pi z\right)$.
$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
\begin{align}
\color{#c00000}{\int_{0}^{1}{\dd x \over \pars{x^{2} - x^{3}}^{1/3}}}&
=\int_{\infty}^{1}{-\dd x/x^{2} \over \pars{1/x^{2} - 1/x^{3}}^{1/3}}
=\int_{1}^{\infty}{\dd x \over x\pars{x - 1}^{1/3}}
=\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}
\end{align}
Use the following contour with
$\ds{z^{-1/3} = \verts{z}^{-1/3}\expo{-\ic\phi\pars{z}/3}\,,\qquad
0 < \phi\pars{z} < 2\pi}$:
\begin{align}
\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}
&=2\pi\ic\expo{-\pi\ic/3}
-\int_{\infty}^{0}{x^{-1/3}\expo{-2\pi\ic/3}\,\dd x \over x + 1}
=2\pi\ic\expo{-\pi\ic/3}
+\expo{-2\pi\ic/3}\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}
\end{align}
\begin{align}
\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}&=
2\pi\ic\,{\expo{-\pi\ic/3} \over 1 - \expo{-2\pi\ic/3}}
=
2\pi\ic\,{1 \over \expo{\pi\ic/3} - \expo{-\pi\ic/3}}
={2\pi\ic \over 2\ic\sin\pars{\pi/3}}={\pi \over \sin\pars{\pi/3}}
\\[3mm]&={\pi \over \root{3}/2} = {2\root{3} \over 3}\,\pi
\end{align}
$$
\color{#66f}{\large\int_{0}^{1}{\dd x \over \pars{x^{2} - x^{3}}^{1/3}}
={2\root{3} \over 3}\,\pi} \approx 3.6276
$$
Best Answer
The function to be integrated (and any related function that will ultimately solve the problem) hates and wants to avoid the points $0,1$, so the pink contour would be a first guess to apply complex analysis (C.A. for short below):
Let $a$ be $a=2\pi i$. We quickly find a polynomial $P=P_a$, so that $P(X+a)-P(X)=X^2$, this is $$ P(X) = \frac 1{3a}\left(x^3-\frac 32ax^2+\frac 12a^2x\right) \ . $$ Then $$ \int_{\color{magenta}{\text{Pink contour}}} \frac {P(\ln z)}{(z-1)^2} =2\pi\; i\sum_{r\text{ Residue}}\operatorname{Res}_{z=r} \frac {P(\ln z)}{(z-1)^2} =0 \ . $$ Now we let $R$ go to $+\infty$.
The upper blue segment produces an integral of the shape $$\displaystyle \int_{0+\varepsilon}^\infty \frac{P(\ln x)}{(x-1)^2}\; dx\ , $$
then the semicircle does not contribute, since for $z=R e^{it}$ the numerator contributes with $\ln^3 R+O(1)$, the denominator is $O(R^{-2})$, the whole fraction is $\sim R^{-2}\ln R$, and the contour is in $O(R)$,
and finally the lower pink segment produces an integral of the shape $$\displaystyle \int_\infty^{0+\varepsilon} \frac{P(\ln x+2\pi i)}{(x-1)^2}\; dx\ , $$ if we succeed to push it beyond the holes.
Here i have to appologize for using $\ln$ for the real function $\log$, and also for the used branch of the logarithm, which is the reason for getting $\ln(x+2\pi i)$ by the monodromy of the logarithm when the contour comes back.
By the choice of $P$, the two integrals come together to build (up to sign) $$ \int_{0+\varepsilon}^\infty \frac{P(\ln x+2\pi i)-P(\ln x)}{(x-1)^2}\; dx = \int_{0+\varepsilon}^\infty \frac{x^2}{(x-1)^2}\; dx \ , $$ and this value comes then during "pushing" exactly from the holes in $0,1$. Let us compute the two contributions. Let $C(z_0,\varepsilon)$ the contour / the circle centered in $z_0\in\Bbb C$ with radius $\varepsilon>0$. Then using the substitution $z=z_0+\varepsilon e^{it}$, $t\in[0,2\pi]$, $$ \left| \int_{C(0,\varepsilon)} \frac {\ln^k z} {(z-1)^2}\; dz \right| \sim 2\pi\epsilon\cdot |\ln\varepsilon+O(1)|^k \in O(\varepsilon^{1/2}) $$ produces no contribution, and $$ \begin{aligned} \int_{C(1,\varepsilon)} \frac {\ln^k z} {(z-1)^2}\; dz &= \int_0^{2\pi} \frac {\displaystyle\ln^k(1+\varepsilon e^{it})} {\displaystyle(1+\varepsilon e^{it}-1)^2} \; i\varepsilon e^{it}\; dt \\ &= \int_0^{2\pi} \frac {\displaystyle \left( \frac 11\varepsilon e^{it} -\frac 12\varepsilon^2 e^{2it} +\frac 13\varepsilon^3 e^{3it} \pm\dots \right)^k} {\displaystyle \varepsilon^2 e^{2it}} \; i\varepsilon e^{it}\; dt \\ &\to \begin{cases} 2\pi i = a&\text{ for }k=1\ ,\\ 0&\text{ for }k\ge 2\ , \end{cases} \qquad \text{ when } \varepsilon\to0\ . \end{aligned} $$ The coefficient of $x$ in $P$ was $\frac a6$, so we get the contribution from $1$ in the form (well, computations were done up to sign) $$ \pm \frac 16a^2=\pm\frac 16(2\pi i)^2=\mp \color{blue}{\frac23\pi^2}\ . $$ (We are choosing the plus sign of course, the integral is built w.r.t. a positive function.)
$\color{red}\square$
Job done, but...
This method applies mot-a-mot to the integrals $$ \int_0^\infty \frac{\ln^4 x}{(x-1)^2}\; dx\ ,\qquad \int_0^\infty \frac{\ln^6 x}{(x-1)^2}\; dx\ ,\qquad \int_0^\infty \frac{\ln^8 x}{(x-1)^2}\; dx\ ,\qquad \dots $$ in the following way:
Numerical experiments.
Numerically, i prefer to compute the half value, $$ \begin{aligned} J(2n) &= \int_1^\infty \frac{\ln^{2n} x}{(x-1)^2}\; dx \qquad\text{ Substitution: } y=\frac{1}{x},\ x=\frac{1}{y},\ dx=-\frac 1{y^2}\; dy \\ &= \int_0^1 \frac{\ln^{2n} y}{(1-y)^2}\; dy =\frac 12\int_0^\infty \frac{\ln^{2n} x}{(x-1)^2}\; dx \ . \end{aligned} $$ Then the expected value of the integral is numerically validated:
Also for some two bigger values: