The simplest Taylor approximation of a function f which has multiple parameters is:

f(x⃗+Δx⃗)≈f(x⃗)+(∇_x f(x⃗))∙Δx⃗+÷{1}{2!}⋅Δx⃗^T⋅H⋅Δx⃗+...

(where H is the Hesse Matrix)

In Electrodynamics, functions of the following form appear for the electrostatic potential:

f(z⃗):=÷{1}{||z⃗||_2}

An example is:

ϕ(x⃗)=k_1⋅∫÷{ρ(y⃗)}{|x⃗-y⃗|}⋅d³y

How do we calculate the value of this function?

ϕ(x⃗)=k_1⋅∫ρ(y⃗)⋅f(x⃗-y⃗)⋅d³y
Δx⃗:=-y⃗
f(x⃗-y⃗)=?

In order to find out, we develop the function into a Taylor Series with development point x⃗.

f(x⃗-y⃗)≈f(x⃗)+(∇_x f(x⃗))∙(-y⃗)+÷{1}{2!}⋅(-y⃗)^T⋅H⋅(-y⃗)+...
f(x⃗-y⃗)≈f(x⃗)+(÷{∂}{∂x^i} f(x⃗))⋅(-y^i)+÷{1}{2!}⋅(-y^i)⋅(÷{∂}{∂x^i} ÷{∂}{∂x^j} f(x⃗))⋅(-y^j)+...

We now actually plug the function f(z⃗):=÷{1}{|z⃗|} in:

f(x⃗-y⃗)≈÷{1}{|x⃗|}+(÷{∂}{∂x^i} ÷{1}{|x⃗|})⋅(-y^i)+÷{1}{2!}⋅(-y^i)⋅(÷{∂}{∂x^i} ÷{∂}{∂x^j} f(x⃗))⋅(-y^j)+...

So, what is ÷{∂}{∂x^i} ÷{1}{|x⃗|}?

÷{∂}{∂x^i} ÷{1}{|x⃗|}=÷{∂}{∂x^i} (x^k⋅x^k)^{-÷{1}{2}}
÷{∂}{∂x^i} ÷{1}{|x⃗|}=-÷{1}{2}⋅{(x^k⋅x^k)}^{-÷{3}{2}}⋅2⋅x^i
÷{∂}{∂x^i} ÷{1}{|x⃗|}=-{(x^k⋅x^k)}^{-÷{3}{2}}⋅x^i

Hence,

f(x⃗-y⃗)≈÷{1}{|x⃗|}+(-{(x^k⋅x^k)}^{-÷{3}{2}})⋅x^i⋅(-y^i)+÷{1}{2!}⋅(-y^i)⋅(÷{∂}{∂x^i} ÷{∂}{∂x^j} f(x⃗))⋅(-y^j)+...
f(x⃗-y⃗)≈÷{1}{|x⃗|}+{(x^k⋅x^k)}^{-÷{3}{2}}⋅x^i⋅y^i+÷{1}{2!}⋅(-y^i)⋅(÷{∂}{∂x^i} ÷{∂}{∂x^j} f(x⃗))⋅(-y^j)+...

Let r:=√{x^k⋅x^k}. Then,

f(x⃗-y⃗)≈÷{1}{|x⃗|}+r^{-3}⋅x^i⋅y^i+÷{1}{2!}⋅(-y^i)⋅(÷{∂}{∂x^i} ÷{∂}{∂x^j} f(x⃗))⋅(-y^j)+...
f(x⃗-y⃗)≈÷{1}{|x⃗|}+÷{x^i⋅y^i}{r^3}+÷{1}{2!}⋅(-y^i)⋅(÷{∂}{∂x^i} ÷{∂}{∂x^j} f(x⃗))⋅(-y^j)+...

So, what is ÷{∂}{∂x^j} ÷{∂}{∂x^i} ÷{1}{|x⃗|}?

÷{∂}{∂x^j} ÷{∂}{∂x^i} ÷{1}{|x⃗|}=÷{∂}{∂x^j} (-{(x^k⋅x^k)}^{-÷{3}{2}}⋅x^i)
÷{∂}{∂x^j} ÷{∂}{∂x^i} ÷{1}{|x⃗|}=(÷{∂}{∂x^j} (-{(x^k⋅x^k)}^{-÷{3}{2}}))⋅x^i+(-{(x^k⋅x^k)}^{-÷{3}{2}})⋅δ^{ij}
(÷{∂}{∂x^j} (-{(x^k⋅x^k)}^{-÷{3}{2}}))=÷{3}{2}⋅{(x^k⋅x^k)}^{-÷{5}{2}}⋅2⋅x^j
(÷{∂}{∂x^j} (-{(x^k⋅x^k)}^{-÷{3}{2}}))=3⋅{(x^k⋅x^k)}^{-÷{5}{2}}⋅x^j
÷{∂}{∂x^j} ÷{∂}{∂x^i} ÷{1}{|x⃗|}=3⋅{(x^k⋅x^k)}^{-÷{5}{2}}⋅x^j⋅x^i+(-{(x^k⋅x^k)}^{-÷{3}{2}})⋅δ^{ij}
÷{∂}{∂x^j} ÷{∂}{∂x^i} ÷{1}{|x⃗|}=3⋅r^{-5}⋅x^j⋅x^i-r^{-3}⋅δ^{ij}
÷{∂}{∂x^j} ÷{∂}{∂x^i} ÷{1}{|x⃗|}=3⋅÷{x^i⋅x^j}{r^5}-÷{δ^{ij}}{r^3}

Hence,

f(x⃗-y⃗)≈÷{1}{|x⃗|}+÷{x^i⋅y^i}{r^3}+÷{1}{2!}⋅(-y^i)⋅(3⋅÷{x^i⋅x^j}{r^5}-÷{δ^{ij}}{r^3})⋅(-y^j)+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i⋅y^i}{r^3}+÷{1}{2!}⋅y^i⋅(3⋅÷{x^i⋅x^j}{r^5}-÷{δ^{ij}}{r^3})⋅y^j+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i⋅y^i}{r^3}+÷{1}{2!}⋅3⋅÷{y^i⋅x^i⋅x^j⋅y^j}{r^5}-÷{y^i⋅δ^{ij}⋅y^j}{r^3}+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i⋅y^i}{r^3}+÷{1}{2!}⋅(3⋅÷{y^i⋅x^i⋅x^j⋅y^j}{r^5}-÷{r^2⋅y^i⋅δ^{ij}⋅y^j}{r^5})+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i⋅y^i}{r^3}+÷{1}{2!}⋅÷{y^i⋅y^j}{r^5}⋅(3⋅x^i⋅x^j-r^2⋅δ^{ij})+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i}{r^3}⋅y^i+÷{1}{2!}⋅÷{1}{r^5}⋅(3⋅y^i⋅y^j⋅x^i⋅x^j-x^k⋅x^k⋅y^i⋅y^j⋅δ^{ij})+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i}{r^3}⋅y^i+÷{1}{2!}⋅÷{1}{r^5}⋅(3⋅y^i⋅y^j⋅x^i⋅x^j-x^i⋅x^j⋅δ^{ij}⋅y^k⋅y^k)+...
f(x⃗-y⃗)≈÷{1}{r}+÷{x^i}{r^3}⋅y^i+÷{1}{2!}⋅÷{x^i⋅x^j}{r^5}⋅(3⋅y^i⋅y^j-δ^{ij}⋅y^k⋅y^k)+...

After plugging this into the original integral:

ϕ(x⃗)=k_1⋅∫÷{ρ(y⃗)}{|x⃗-y⃗|}⋅d³y
ϕ(x⃗)=k_1⋅(÷{1}{r}⋅Q_{ges}+÷{x^i}{r^3}⋅p^i+÷{x^i⋅x^j}{r^5}⋅Q^{ij}+...)

With:

Q_{ges}:=∫ρ(y⃗)⋅d³y Monopole momentum
p^i:=∫ρ(y⃗)⋅y^i⋅d³y Dipole momentum
Q^{ij}:=÷{1}{2!}⋅∫ρ(y⃗)⋅(3⋅y^i⋅y^j-y^k⋅y^k⋅δ^{ij})⋅d³y Quadrupole momentum