Solution (source code)

= Solution

Let $\theta$ be the angle between $r$ and the positive $z$-axis. By superposition of the <electric potential> of point charges,
$$
\Phi(r)=\frac{Q}{4\pi\epsilon_0}
\left(\frac N r-\frac1{|r-d\widehat z|}
-\frac M{|r+d\widehat z|}\right).
$$
For $r>d$, the <Legendre polynomial> expansion gives
$$
\frac1{|r\mp d\widehat z|}
=\frac1r\left[1\pm\frac d r\cos\theta
+\frac{d^2}{2r^2}(3\cos^2\theta-1)+O(r^{-3})\right].
$$
Therefore
$$
\boxed{\Phi(r,\theta)=\frac{Q}{4\pi\epsilon_0}\left[
\frac{N-M-1}{r}
+\frac{(M-1)d\cos\theta}{r^2}
-\frac{(M+1)d^2(3\cos^2\theta-1)}{2r^3}
+O(r^{-4})\right]}.
$$
The three displayed terms are respectively the monopole, dipole, and quadrupole terms of the <electric multipole expansion>. In particular, the total charge is $Q(N-M-1)$ and the <electric dipole moment> is $Qd(M-1)\widehat z$. This is the <multipole expansion of three collinear charges>.