I think you should take a look at your derivative again to make sure it is right (it is not quite, but close). You are on the right track, you have the form $\frac{1}{1+y}$, you want $\frac{1}{1-y}$, what do you have to do to your $y$ term to get this? When someone told you to integrate, they mean integrate the power series. You know you can do this in the interval of convergence, right? You are right, once you integrate it you will get $\arctan$ again, but you will also have the integral of the power series on one side which is then the power series for $\arctan$!
Your insistence "let's keep everything in real variables, please" is precisely the problem: the cause of the finite radius of convergence is due to the function's behavior in $\mathbf C$, not $\mathbf R$.
A much simpler example than $\arctan x$ is $1/(1+x^2)$, which is defined and infinitely differentiable on the whole real line but its power series at $0$ (a geometric series with $-x^2$ in place of $x$) has radius of convergence $1$, not $\infty$. To use your language, "there is an obvious discontinuity that you bump into as you work your way out from the center," namely at $x = \pm i$ where the function blows up. In fact, if you expand $1/(1+x^2)$ into a power series at a real number $a$, not necessarily at $0$, the radius of convergence will be $\sqrt{a^2+1} = |a-i|$ -- the distance from the center out to $i$. This phenomenon is bewildering if you refuse to use complex numbers and extremely clear if you use them. Choose wisely.
If $f(x)$ is a rational function in reduced form with a nonconstant denominator and its denominator does not vanish at $a$, its power series at $a$ has radius of convergence $|a-\rho|$ where $\rho$ is the root of the denominator in $\mathbf C$ that is closest to $a$. This simple geometric result can not be explained in terms of real variables if the roots of the denominator are not all real.
To reinforce how poorly the real numbers are compared to the complex numbers as a predictive tool for the radius of convergence, there are functions $\mathbf R \rightarrow \mathbf R$ that are infinitely differentiable on the whole real line but their power series at each real number $a$ has radius of convergence zero for all $a$ in $\mathbf R$.
Strictly speaking, the real numbers have enough information in principle to compute the radius of convergence $R$ of a power series $\sum c_n(x-a)^n$ with all real coefficients using Hadamard's formula $1/R = \varlimsup\limits_{n\to\infty} \sqrt[n]{|c_n|}$, but this formula is often not feasible to compute in practice.
Best Answer
Yes, there is indeed a solution by differentiation + integration:
Differentiate your expression:
$$2x \arctan(x) +1$$
replace $\arctan(x)$ by its expansion, then integrate it term by term... taking into account the integration constant in such a way that there is an agreement in $x=0$ of the initial expression and your expansion.