Let $X_i, X_j$ be two random variables that can each assume the values $\zeta_1, ..., \zeta_m$. Then my book claims $E(X_iX_j) = \sum_{k = 1}^m\sum_{l = 1}^m\zeta_k\zeta_lP(X_i = \zeta_k \textrm{ and } X_j = \zeta_l$).
I don't get how they acquired this result. To my understanding $E(X_iX_j) = \sum_{k = 1}^m\sum_{l = 1}^m\zeta_k\zeta_lP(X_i X_j = \zeta_k\zeta_l)$. So why is $P(X_iX_j = \zeta_k\zeta_l) = (X_i = \zeta_k \textrm{ and } X_j = \zeta_l)$?
We do have: $$\mathsf E(X_iX_j)=\sum_{z\in Z} z\mathsf P(X_iX_j=z)\tag1$$ where $Z:=\{\zeta_k\zeta_l\mid k,l\in\{1,\dots,m\}\}$
Every $z\in Z$ can be written as $z=\zeta_k\zeta_l$ for some pair $k,l\in\{1,\dots,m\}$ but this value must be encountered only once in the RHS of $(1)$.
According to your (wrong) understanding it can occur more than once.
The claim of your book is correct. There we also meet values $z=\zeta_k\zeta_l\in Z$ and they can occur more than once. However they are not multiplied by $\mathsf P(X_iX_j=z)$ but by $\mathsf P(X_i=\zeta_k\wedge X_j=\zeta_l)$.
Observe that for $z\in Z$: $$P(X_iX_j=z)=\sum_{\zeta_k\zeta_l=z} P(X_i=\zeta_k,X_j=\zeta_l)\tag2$$
If the RHS of $(2)$ is substituted in $(1)$ then we eventually will arrive at:$$\mathsf EX_iX_j=\sum_{k=1}^m\sum_{l=1}^m\zeta_k\zeta_l\mathsf P(X_i=\zeta_k\wedge X_j=\zeta_l)\tag3$$i.e. the claim of your book.
Your second formula is "double counting". Take for instance $X_1=B_1+1$ and $X_2=B_2+1$, with $B_1$ and $B_2$ both Bernoulli distributed with the same probability $p$ and independent from each other. Then
$$E(X_1X_2) = 1\cdot 1 \cdot P(X_1=1 \text{ and } X_2=1) + 1\cdot 2 P(X_1=1 \text{ and } X_2=2) \\ + 2\cdot 1 P(X_1=2 \text{ and } X_2=1) + 2\cdot 2 P(X_1=2 \text{ and } X_2=2) = (1-p)^2 + 4p(1-p)+4p^2$$
on the other hand, if your second formula were true you would have
$$E^*(X_1X_2) = 1\cdot 1 \cdot P(X_1 X_2=1) + 1\cdot 2 P(X_1 X_2=2) \\ + 2\cdot 1 P(X_1X_2=2) + 2\cdot 2 P(X_1 X_2=4) = (1-p)^2 + 8p(1-p)+4p^2$$
The middle terms have been double counted because $P(X_1 X_2=2)=P(X_1=1 \text{ and } X_2=2)+P(X_1=2 \text{ and } X_2=1)$. So the second formula is incorrect.
This stuff always confuses me, so I go back to the basic derivation. The relevant theorem is:
Let $X=(X_{1},...,X_{n})$ be an n-dimensional discrete random variable with Range $R_X$ and density $f_X$, and let $g:\mathbb{R}^n \longrightarrow \mathbb{R}$ be a measurable function. Then: $$\mathbb{E}(g(X)) = \sum_{r\in R_X} g(r)f_X(r)$$
In this case we have the 2-dimensional discrete random variable $(X,Y)$. Let $g(x,y)=xy$, then from the theorem we have: $$\mathbb{E}(XY) = \mathbb{E}(g(X,Y)) = \sum_{(x,y)\in R_{X,Y}} g(x,y)f_{X,Y}(x,y)=\\\sum_{x\in R_X}\sum_{y\in R_Y} xy\cdot\mathbb{P}(X=x,Y=y)$$
