Given two vectors $\mathbf x\in\mathbb{R}^3$ and $\mathbf x'\in\mathbb{R}^3$. Assume: $x = |\mathbf x|$ and $x' = |\mathbf x'|$.
Show that: $$ \frac{1}{|\mathbf x - \mathbf x'|} = \frac{1}{\sqrt{x^2 + x'^2 - 2xx'\cos(\mathbf x, \mathbf x')}} = \sum_{i=0}^\infty\frac{x^i_<}{x^{i+1}_>} P_i(\cos(\mathbf x, \mathbf x')) $$
Where $x_>$ and $x_<$ are defined: $$ x_> = \begin{cases} x & \text{if } x > x' \\ x' & \text{if } x < x' \end{cases} \quad\quad x_< = \begin{cases} x & \text{if } x < x' \\ x' & \text{if } x > x' \end{cases} $$
and where $P_i(x)$ is the Legendre polynomials. A representation: $$ P_i(x) = \frac{1}{2^i i!}\frac{d^i}{dx^i}\left(x^2-1\right)^i $$
I think this is done using Taylor expansion. I tried but I failed. Actually, I have no clue where to begin to prove this.
Taylor expansion won't get you where you wish to go. Start with the general solution for Laplace's equation, which involves spherical harmonics: $\Delta \Phi = 0$ solutions are $$\Phi(r,\theta,\phi) = \sum_{l=0}^\infty\sum_{m=-l}^l [A_{lm}r^l+B_{lm}r^{-(l+1)}]N_{lm}P_{lm}(cos\theta)e^{im\phi}$$ Where $N_{lm}$ is a normalizing factor which may be absorbed into $A_{lm}$ and $B_{lm}$.
As $\frac{1}{|\mathbf{r}-\mathbf{r^{'}}|}$ is a solution to the Laplace equation, it can be written as above. Obviously then, only $m=0$ applies. So it remains to determine $A_l$ and $B_l$ in
$$\frac{1}{|\mathbf{r}-\mathbf{r^{'}}|}=\sum_{l=0}^\infty[A_lr^l+B_lr^{-(l+1)}]P_{l}(cos\theta)$$
First let's look at an $\mathbf{r^{'}}$ vector along the z-axis, and also first take $\mathbf{r}$ aligned with $\mathbf{r^{'}}$, so that $\theta=0$, and (remember that $P_{l}(1)=1$ for all $l$)
$$\frac{1}{|r-r^{'}|}=\sum_{l=0}^\infty[A_lr^l+B_lr^{-(l+1)}]$$
It is also clear that
$$\frac{1}{|r-r^{'}|}=\frac{1}{r-r^{'}}=\frac{1}{r}\frac{1}{1-\frac{r^{'}}{r}}=\frac{1}{r}\sum_{l=0}^\infty(\frac{r^{'}}{r})^l=\sum_{l=0}^\infty(r^{'})^l\cdot r^{-(l+1)}$$ when $r>r^{'}$
and
$$\frac{1}{|r-r^{'}|}=\frac{1}{r^{'}-r}=\frac{1}{r^{'}}\frac{1}{1-\frac{r}{r^{'}}}=\frac{1}{r^{'}}\sum_{l=0}^\infty(\frac{r}{r^{'}})^l=\sum_{l=0}^\infty({r^{'}})^{-(l+1)}\cdot r^l$$ when $r<r^{'}$
So the solution splits in two functions: one for $r>r^{'}$ where $A_l=0$ and $B_l=(r^{'})^l$, and one for $r<r^{'}$ where $A_l=(r^{'})^{-(l+1)}$ and $B_l=0$.
Clearly, then, for other $\mathbf{r}$ that are not aligned with $\mathbf{r^{'}}$, the values of the factors $A_l$ and $B_l$ are the same but the series expansion now includes factors $P_{l}(cos\theta)$, where $\theta$ is the angle between $\mathbf{r}$ and $\mathbf{r^{'}}$ (which obviously is the same as the actual $\theta$, the angle between $\mathbf{r}$ and the z-axis).
Moving to $\mathbf{r^{'}}$ not along the z-axis, the same expansion will hold, as rotating both $\mathbf{r^{'}}$ and $\mathbf{r}$ in the same fashion does not change the angle between them, nor their length.
Hence
$$\frac{1}{|\mathbf{r}-\mathbf{r^{'}}|}=\sum_{l=0}^\infty\frac{r_{<}^l}{r_{>}^{l+1}}P_{l}(cos\theta)$$