Using power series to solve differential equations of high order I have to solve this different equation with the help of power series
$y''+y'+y = x+\frac{1}{3}x^{3}+\frac{1}{5}x^{5}...$
Let $y = \sum_{0}^{\infty}a_{n}x^{n}$
$y' = \sum_{1}^{\infty}n\cdot a_{n}x^{n-1}$
$y'' = \sum_{2}^{\infty}n\cdot(n-1)\cdot a_{n}x^{n-2}$
Now the exponents in $x^{n}$ of all the derivatives must be the same for $y'$ $n\rightarrow n+1$
$y' = \sum_{0}^{\infty}(n+1)\cdot a_{n+1}x^{n}$
and for $y''$ $n\rightarrow n+2$
$y'' = \sum_{2}^{\infty}(n+2)\cdot(n+1)\cdot a_{n+2}x^{n}$ $\rightarrow \sum_{0}^{\infty}[a_{n}+(n+1)a_{n+1}+(n+2)(n+1)a_{n+2}]x^n =$
Now if we had 0 instead of $x+\frac{1}{3}x^{3}+\frac{1}{5}x^{5}...$ then we could find the first 4 non zero coefficients but now how do I continue?

The homogeneous equation
$$u''+u'+u=0$$
Is solved very easily with exponentials, so let's find a particular solution to
$$u''(x)+u'(x)+u(x)=\operatorname{arctanh}(x)=\sum_{k=1}^\infty \frac{z^{2k-1}}{2k-1}$$
Using a power series. First, we remark that
$$\sum_{k=1}^\infty \frac{1}{2k-1}x^{2k-1}=\sum_{n=1}^\infty \frac{1-(-1)^n}{2n}x^n$$
Now,
$$u(x)=\sum_{n=0}^\infty a_nx^n  \\ u'(x)=\sum_{n=0}^\infty (n+1)a_{n+1}x^n \\ u''(x)=\sum_{n=0}^\infty (n+1)(n+2)a_{n+2}x^n$$
$$\sum_{n=0}^\infty \big(a_n+(n+1)a_{n+1}+(n+1)(n+2)a_{n+2}\big)x^n=\sum_{n=1}^\infty \frac{1-(-1)^n}{2n}x^n\tag{1}$$
Because the right hand side has no constant term, the left must be the same:
$$a_0+a_{1}+2a_2=0\tag{2}$$
We can give this whatever initial values we want, so let's pick for instance $u'(0)=u(0)=0$.
From this we get
$$a_0=a_1=0$$
Which, from $(2)$, gives us $a_2=0$ as well.
Going back to $(1)$ we as well we get the recurrence
$$a_n+(n+1)a_{n+1}+(n+1)(n+2)a_{n+2}=\frac{1-(-1)^n}{2n}$$
Rearranging,
$$(n+1)(n+2)a_{n+2}=\frac{1-(-1)^n}{2n}-(n+1)a_{n+1}-a_n \\ a_{n+2}=\frac{1-(-1)^n}{2n(n+1)(n+2)}-\frac{1}{n+2}a_{n+1}-\frac{1}{(n+1)(n+2)}a_n $$
Shifting the index,
$$a_n=\frac{1-(-1)^n}{2n(n-1)(n-2)}-\frac{1}{n}a_{n-1}-\frac{1}{n(n-1)}a_{n-2}$$
From this we can easily generate the coefficients quickly using python, for instance
import numpy as np
nmax=10
a=np.zeros(nmax)
a[0]=0
a[1]=0
a[2]=0
for n in range(3,nmax):
    a[n]=(1-(-1)**n)/(2*n*(n-1)*(n-2))-a[n-1]/n-a[n-2]/(n*(n-1))
print(str(a))

Returns the output
[ 0.          0.          0.          0.16666667 -0.04166667  0.01666667
 -0.00138889  0.00456349 -0.00054563  0.00198137]

These are the values of $(a_n)_{0\leq n\leq 9}$.
Then you can plot it,
def f(x):
    return sum(a[n]*x**n for n in range(0,nmax))
x=np.arange(-0.99,1,0.010)
plt.plot(x,f(x))
plt.xlim([-1,1])
plt.ylim([-0.25,0.25])
plt.title("Approximate solution of $u''(x)+u'(x)+u(x)=tanh^{-1}(x)$ with $u(0)=u'(0)=0$ up to $O(x^{10})$")

Output:

