Maclaurin series for $\arctan^{2}(x)$

Using the Cauchy Product Formula, $$ \begin{align} \frac{\mathrm{d}}{\mathrm{d}x}\arctan^2(x) &=2\frac{\arctan(x)}{1+x^2}\\ &=2\left(x-\frac{x^3}3+\frac{x^5}5-\frac{x^7}7+\cdots\right)\left(1-x^2+x^4-x^6+\cdots\right)\\ &=2\left((1)x-\left(1+\tfrac13\right)x^3+\left(1+\tfrac13+\tfrac15\right)x^5-\left(1+\tfrac13+\tfrac15+\tfrac17\right)x^7+\cdots\right)\\ &=2\sum_{k=1}^\infty(-1)^{k-1}\left(H_{2k}-\tfrac12H_k\right)x^{2k-1} \end{align} $$ Therefore, $$ \arctan^2(x)=\sum_{k=1}^\infty\frac{(-1)^{k-1}}k\left(H_{2k}-\tfrac12H_k\right)x^{2k} $$

We can try to obtain the series in the following way:

$$f(x)=\arctan^2 x=x^2 \int_0^1 \int_0^1 \frac{du~dv}{(1+x^2u^2)(1+x^2v^2)}$$

It's easier to consider:

$$g(x)=\int_0^1 \int_0^1 \frac{du~dv}{(1+x^2u^2)(1+x^2v^2)}$$

Let's use partial fractions:

$$\frac{1}{(1+x^2u^2)(1+x^2v^2)}=\frac{u^2}{(u^2-v^2)(1+x^2u^2)}-\frac{v^2}{(u^2-v^2)(1+x^2v^2)}$$

We obtain a sum of two singular integrals, which however, can both be formally expanded into a series:

$$g(x)=\int_0^1 \int_0^1 \left(\frac{u^2}{(u^2-v^2)(1+x^2u^2)}-\frac{v^2}{(u^2-v^2)(1+x^2v^2)} \right) du ~dv=$$

$$g(x)=\sum_{n=0}^\infty (-1)^n x^{2n} \int_0^1 \int_0^1 \left(\frac{u^{2n+2}}{u^2-v^2}-\frac{v^{2n+2}}{u^2-v^2} \right) du ~dv$$

Now:

$$g(x)=\sum_{n=0}^\infty (-1)^n x^{2n} \int_0^1 \int_0^1 \frac{u^{2n+2}-v^{2n+2}}{u^2-v^2} du ~dv$$

Obviously, every integral is finite now, and we can write:

$$\int_0^1 \int_0^1 \frac{u^{2n+2}-v^{2n+2}}{u^2-v^2} dv ~du=2 \int_0^1 \int_0^u \frac{u^{2n+2}-v^{2n+2}}{u^2-v^2} dv ~du= \\ = 2 \sum_{k=0}^\infty \int_0^1 \int_0^u u^{2n} \left(1-\frac{v^{2n+2}}{u^{2n+2}} \right) \frac{v^{2k}}{u^{2k}} dv ~du =2 \sum_{k=0}^\infty \int_0^1 \int_0^1 u^{2n+1} \left(1-t^{2n+2} \right) t^{2k} dt ~du= \\ = 2 \sum_{k=0}^\infty \int_0^1 \left(\frac{1}{2k+1}-\frac{1}{2k+2n+3} \right) u^{2n+1}~du= 2\sum_{k=0}^\infty \frac{1}{(2k+1)(2k+2n+3)}$$

So we get:

$$g(x)=\frac{1}{2} \sum_{n=0}^\infty (-1)^n x^{2n} \sum_{k=0}^\infty \frac{1}{(k+\frac{1}{2})(k+n+\frac{3}{2})}$$

The inner series converges for all $n$, and we can formally represent it as a difference of two divergent harmonic series, which, after some manipulations, should give us the same result as robjohn obtained.

I suppose, a kind of closed form for the general term can also be given in terms of digamma function:

$$g(x)=\frac{1}{2} \sum_{n=0}^\infty (-1)^n \frac{\psi \left(n+\frac32 \right)-\psi \left(\frac12 \right)}{n+1} x^{2n}$$

Which makes:

$$\arctan^2 x=\frac{x^2}{2} \sum_{n=0}^\infty (-1)^n \frac{\psi \left(n+\frac32 \right)-\psi \left(\frac12 \right)}{n+1} x^{2n}$$

Which is essentially the same as robjohn's answer.

For the sake of an alternative method, here's a double series.

For $|x|<1$, we have $$\arctan x=\sum_{k\geq0}\frac{(-1)^kx^{1+2k}}{1+2k}$$ As you noted, $$\arctan^2x=2\int\frac{\arctan x}{x^2+1}\mathrm{d}x$$ So, assuming $|x|<1$, $$\arctan^2x=2\int\frac1{x^2+1}\sum_{k\geq0}\frac{(-1)^kx^{1+2k}}{1+2k}\mathrm{d}x$$ $$\arctan^2x=2\sum_{k\geq0}\frac{(-1)^k}{1+2k}\int\frac{x^{2k}}{1+x^2}x\mathrm{d}x$$ Now we focus on $$I=\int\frac{x^{2k+1}}{1+x^2}\mathrm{d}x$$ We recall that for $|x|<1$, $$\frac1{1+x^2}=\frac12\sum_{n\geq0}i^n\big(1+(-1)^n\big)x^n$$ Hence $$I=\frac12\sum_{n\geq0}i^n\big(1+(-1)^n\big)\int x^{2k+n+1}\mathrm{d}x$$ $$I=\frac12\sum_{n\geq0}\frac{i^n\big(1+(-1)^n\big)}{2k+n+2}x^{2k+n+2}$$ Then we have our (pretty inefficient) result: $$\arctan^2x=\sum_{k\geq0}\frac{(-1)^k}{1+2k}\sum_{n\geq0}\frac{i^n\big(1+(-1)^n\big)}{2k+n+2}x^{2k+n+2}$$ $$\arctan^2x=\sum_{n,k\in\Bbb N}\frac{(-1)^{k+\frac{n}2}\big(1+(-1)^n\big)}{(1+2k)(2+2k+n)}x^{2k+n+2}$$