Find $f(x)$ such that: $ f'(x) + f(x^2) = 2x + 1 $ How to find the function $f(x)$ that is derivable on $\mathbb{R}$ and satisfies the equation:
$$ f'(x) + f(x^2) = 2x + 1 \text{ } \text{ } \forall x \in \mathbb{R}$$

My attempt:
Substitute $-x$ for $x$ in the equation, we have a system of 2 equations :
  \begin{cases} 
f'(x) + f(x^2) &= 2x + 1\\ 
f'(-x) + f(x^2) &= -2x + 1
\end{cases} 
  Take the difference of them, we have: 
  $$ f'(x) - f'(-x) =4x $$ 
  Integrate both sides: 
  $$ f(x) + f(x) = 2x^2 + C_1 $$
  Therefore, we have:
  $$ f(x) = x^2 + C_2 $$

However, this method results in an invalid solution to the original equation. I wonder whether there is another way to solve this problem or why my solution is wrong.
Thanks in advance.
 A: We can write a power series for $f(x)$ as $$f(x) = \sum_{n=0}^\infty a_nx^n$$
Then we get that $$f'(x) + f(x^2) = \sum_{n=0}^{\infty}(n+1)a_{n+1}x^{n} + \sum_{n=0}^{\infty} a_nx^{2n} = 2x+1$$
Breaking this into odd and even powers of $x$, this can be rewritten as $$\sum_{n=0}^{\infty}((2n+1)a_{2n+1} + a_n)x^{2n} + \sum_{n=0}^{\infty}(2n+2)a_{2n+2}x^{2n+1} = 2x+1$$
From the second sum, we can get that $(2n+2)a_{2n+2} = 0 \to a_{2n+2} = 0$ for $n \ge 1$ and $a_2 = 1$. Then from the first sum, we have these relations. 
$$a_1 + a_0 = 1$$
$$3a_3 + a_1 = 0$$
$$5a_5 + a_2 = 0$$
$$7a_7 + a_3 = 0$$
$$9a_9 + a_4 = 0$$
$$...$$
Because $a_4, a_6, a_8, a_{10}, ...$ are $0$, it is also true that $a_9, a_{13}, a_{17}, a_{21}, ..$ are also $0$. Because those are $0$, it is also true that $a_{19}, a_{27}, a_{35}, ...$ are also $0$. By induction, the only ones that are not $0$ are the powers of $x$ where the exponent is either of the form $2^n-1, n \ge 0$ or $3 \cdot 2^n - 1, n \ge 0$.
From $a_2 = 1$, we can get that $a_5 = -1/5$, which leads to $a_{11} = 1/55$. In general, $$a_{3 \cdot 2^n - 1} = \frac{(-1)^n}{\prod_{k=1}^{n} (3 \cdot 2^k - 1)}$$
Similarly, given $a_1$, we can get that $a_0 = 1-a_1$, and $$a_{2^n-1} = \frac{(-1)^{n+1}}{\prod_{k=1}^n (2^k - 1)} a_1 $$
From this, we have that $$f(x) = (1-a_1) + x^2 + \sum_{n=1}^{\infty} \frac{(-1)^n}{\prod_{k=1}^{n} (3 \cdot 2^k - 1)} x^{{3 \cdot 2^n - 1}} + a_1\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{\prod_{k=1}^n (2^k - 1)} x^{2^n-1}$$
$a_1$ can be anything. Note that the required functional equation only holds true for $-1 < x < 1$ for my $f(x)$, as it only converges for $-1 < x < 1$.
A: Follows a MATHEMATICA script which calculates the series approximation for this ODE. It is convergent for $|x| < 1$
n = 50;
f[x_, n_] := Sum[Subscript[a, k] x^k, {k, 0, n}]
A = Table[Subscript[a, k], {k, 0, n}];
dif = D[f[x, n], x] + f[x^2, n] - 2 x - 1;
equs = Take[CoefficientList[dif, x], {1, n + 1}];
B = equs /. Thread[A -> 0];
M = Grad[equs, A];
solA = LinearSolve[M, -B];

f[x_] := f[x, n] /. Thread[A -> solA]

f'[x] + f[x^2] - 2 x - 1

Plot[f[x], {x, -1, 1}, PlotStyle -> {Thick, Blue}]


