Today I was playing with this function $$ \sqrt{1-x} $$ and I found that if I try to approximate (near $x=0$) such a function with a polynomial of some order then I can find the same coefficients that WolframAlpha find.
Starting with $$ \sqrt{1-x}\sim a_0+a_1x+a_2x^2+a_3x^3+\cdots $$ where $\cdots$ are terms of higher power in $x$, I squared both sides obtaining $$ 1-x\sim a_0^2+a_1^2x^2+a_2^2x^4+a_3^2x^6+2a_0a_1x+2a_0a_2x^2+2a_0a_3x^3+2a_1a_2x^3+2a_1a_3x^4+2a_2a_3x^5+\dots $$ or $$ 1-x\sim a_0^2+2a_0a_1x+(a_1^2+2a_0a_2)x^2+(2a_0a_3+2a_1a_2)x^3+\dots $$ where I choose to ignore terms $x^n$ with $n>3$. In this way I obtain an "equivalence" between two polynomials and imposing that each term of those polynomials are equal one find $$ a_0^2=1,\,2a_0a_1=-1,\,a_1^2+2a_0a_2=0,\,2a_0a_3+2a_1a_2=0. $$ Choosing $a_0=1$ one get $$ \sqrt{1-x}\sim1-\frac{1}{2}x-\frac{1}{8}x^2-\frac{1}{16}x^3+\cdots. $$ Is this method correct? I tried and it works also for $\frac{1}{1-x}$ or others irrational and rational functions! I also would love if someone give some references.
This is not an accident. It is a general way to find the square root of a power series.
Writing it out:
$\begin{array}\\ \sum_{n=0}^{\infty} c_n x^n &=(\sum_{i=0}^{\infty} a_i x^i)^2\\ &=(\sum_{i=0}^{\infty} a_i x^i)(\sum_{j=0}^{\infty} a_j x^j)\\ &=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty} a_ia_j x^{i+j}\\ &=\sum_{n=0}^{\infty}x^n\sum_{i=0}^{n} a_ia_{n-i}\\ \end{array} $
so $c_n =\sum_{i=0}^{n} a_ia_{n-i} $.
For $n=0$, $c_0 = a_0^2$, so $a_0 = \sqrt{c_0}$.
For $n > 0$, if $n = 2m+1$,
$\begin{array}\\ c_n &=c_{2m+1}\\ &=\sum_{i=0}^{2m+1} a_ia_{2m+1-i}\\ &=\sum_{i=0}^{m} a_ia_{2m+1-i}+\sum_{i=m+1}^{2m+1} a_ia_{2m+1-i}\\ &=\sum_{i=0}^{m} a_ia_{2m+1-i}+\sum_{i=0}^{m} a_{i+m+1}a_{2m+1-(i+m+1)}\\ &=\sum_{i=0}^{m} a_ia_{2m+1-i}+\sum_{i=0}^{m} a_{i+m+1}a_{m-i}\\ &=\sum_{i=0}^{m} a_ia_{2m+1-i}+\sum_{i=0}^{m} a_{(m-i)+m+1}a_{i}\\ &=\sum_{i=0}^{m} a_ia_{2m+1-i}+\sum_{i=0}^{m} a_{2m+1-i}a_{i}\\ &=2\sum_{i=0}^{m} a_ia_{2m+1-i}\\ &=2a_0a_{2m+1}+2\sum_{i=1}^{m} a_ia_{2m+1-i}\\ \end{array} $
so $a_{2m+1} =\dfrac{c_{2m+1}-2\sum_{i=1}^{m} a_ia_{2m+1-i}}{2a_0} $.
If $n = 2m$,
$\begin{array}\\ c_n &=c_{2m}\\ &=\sum_{i=0}^{2m} a_ia_{2m-i}\\ &=\sum_{i=0}^{m-1} a_ia_{2m-i}+a_m^2+\sum_{i=m+1}^{2m} a_ia_{2m-i}\\ &=\sum_{i=0}^{m-1} a_ia_{2m-i}+a_m^2+\sum_{i=0}^{m-1} a_{i+m+1}a_{2m-(i+m+1)}\\ &=\sum_{i=0}^{m-1} a_ia_{2m-i}+a_m^2+\sum_{i=0}^{m-1} a_{i+m+1}a_{m-i-1}\\ &=\sum_{i=0}^{m-1} a_ia_{2m-i}+a_m^2+\sum_{i=0}^{m-1} a_{(m-1-i)+m+1}a_{m-(m-1-i)-1}\\ &=\sum_{i=0}^{m-1} a_ia_{2m-i}+a_m^2+\sum_{i=0}^{m-1} a_{2m-i}a_{i}\\ &=2\sum_{i=0}^{m-1} a_ia_{2m-i}+a_m^2\\ &=2a_0a_{2m}+2\sum_{i=1}^{m-1} a_ia_{2m-i}+a_m^2\\ \end{array} $
so $a_{2m} =\dfrac{c_{2m}-a_m^2-2\sum_{i=1}^{m-1} a_ia_{2m-i}}{2a_0} $.
Note that we must have $c_0 \ne 0$ and there are two possible series for each sign in $a_0 = \pm \sqrt{c_0}$.
If $c_0 < 0$, the result requires complex numbers.
If $c_0 = 1$, then $a_0 = \pm 1$.