MAIN OBJECTIVE : To find the electric field due to a cubic body at free space points (except few) and then ensure that the divergence of its electric field in free space is zero
GIVEN : Charge density is unity at all points inside the cubic body AND the vertices of cube are $(-20,-20,-20)$,$(-20,-20,20)$,$(-20,20,-20)$,$(-20,20,20)$,$(20,-20,-20)$,$(20,-20,20)$,$(20,20,-20)$,$(20,20,20)$
For any $(a,b,c)$ not on cube such that $a \neq 20, a \neq -20, b \neq 20, b \neq -20, c \neq 20, c \neq -20$
Part I : Finding $x$-component of electric field
\begin{align} E_{x}(a,b,c)&=\int_{cube} \dfrac{a-x}{[(a-x)^2+(b-y)^2+(c-z)^2]^{3/2}}\ dV\\ &=\int^{20}_{-20} \int^{20}_{-20} \int^{20}_{-20} \dfrac{a-x}{[(a-x)^2+(b-y)^2+(c-z)^2]^{3/2}}\ dx\ dy\ dz\\ &=\int^{20}_{-20} \int^{20}_{-20} \left[ \dfrac{1}{\sqrt{(a-x)^2+(b-y)^2+(c-z)^2}} \right]^{20}_{-20} dy\ dz\\ &=\int^{20}_{-20} \int^{20}_{-20} \left[ \dfrac{1}{\sqrt{(a-20)^2+(b-y)^2+(c-z)^2}}-\dfrac{1}{\sqrt{(a+20)^2+(b-y)^2+(c-z)^2}} \right] dy\ dz\\ &=\int^{20}_{-20}\left[ \ln \left| \sqrt{(a-20)^2+(b-y)^2+(c-z)^2}+y-b\right| - \ln \left|\sqrt{(a+20)^2+(b-y)^2+(c-z)^2}+y-b \right| \right]^{20}_{-20} dz\\ &=\int^{20}_{-20} \left[ \ln \left| \sqrt{(a-20)^2+(b-20)^2+(c-z)^2}+20-b\right| - \ln \left|\sqrt{(a+20)^2+(b-20)^2+(c-z)^2}+20-b \right|\\ - \ln \left|\sqrt{(a-20)^2+(b+20)^2+(c-z)^2}-20-b \right| + \ln \left|\sqrt{(a+20)^2+(b+20)^2+(c-z)^2}-20-b \right| \right]\ dz\\ &=\int^{20}_{-20} \left[ \ln \left| \sqrt{(a-20)^2+(b-20)^2+(c-z)^2}-(b-20)\right| - \ln \left|\sqrt{(a+20)^2+(b-20)^2+(c-z)^2}-(b-20) \right|\\ - \ln \left|\sqrt{(a-20)^2+(b+20)^2+(c-z)^2}-(b+20) \right| + \ln \left|\sqrt{(a+20)^2+(b+20)^2+(c-z)^2}-(b+20) \right| \right]\ dz\\ \end{align}
We see that there are four terms in the integrand of integral w.r.t $z$. First we consider only the first term. According to this integral calculator, its integral w.r.t $z$ along with limits is:
\begin{align} \text{Term 1} &=\left[ (b-20) \left( \ln{\left| \sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - z)^{2} } - \sqrt{(a - 20)^{2} + (b - 20)^{2} } + (c-z) \right|}\\ -\ln{\left| \sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - z)^{2} } - \sqrt{(a - 20)^{2} + (b - 20)^{2} } - (c-z) \right|} \right)\\ -(c-z) \ln\left|\sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - z)^{2}} - (b - 20) \right|\\ + 2 (a - 20) \tan^{-1} \left( \frac{a-20}{c-z}\ .\ \frac{\sqrt{(a - 20)^{2} + (b - 20)^{2}} - \sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - z)^{2}}}{\sqrt{(a - 20)^{2} + (b - 20)^{2}} - (b-20)} \right)\\ - z \right]^{20}_{-20} \end{align}
EDITED FROM HERE
The last term $-z$ gets cancelled off when applying limits. Therefore the integral w.r.t $z$ of Term $1$ computed for upper limit $z=20$ can be written as the following (call equation $1$):
\begin{align} E_{x(1,2,3)} &=(b-20) \left( \ln{\left| \sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - 20)^{2} } - \sqrt{(a - 20)^{2} + (b - 20)^{2} } + (c-20) \right|}\\ -\ln{\left| \sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - 20)^{2} } - \sqrt{(a - 20)^{2} + (b - 20)^{2} } - (c-20) \right|} \right)\\ &-(c-20) \ln\left|\sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - 20)^{2}} - (b - 20) \right|\\ &+ 2 (a - 20) \tan^{-1} \left( \frac{a-20}{c-20}\ .\ \frac{\sqrt{(a - 20)^{2} + (b - 20)^{2}} - \sqrt{(a - 20)^{2} + (b - 20)^{2} + (c - 20)^{2}}}{\sqrt{(a - 20)^{2} + (b - 20)^{2}} - (b-20)} \right) \tag1 \end{align}
By replacing:
$a-20$ with $a+20$ in equation $(1)$
and/or
$b-20$ with $b+20$ in equation $(1)$
and/or
$c-20$ with $c+20$ in equation $(1)$,
we see that there are $7$ more variations of $a \pm 20$,$b \pm 20$,$c \pm 20$
| a | b | c |
|---|---|---|
| a-20 | b-20 | c-20 |
| a-20 | b-20 | c+20 |
| a-20 | b+20 | c-20 |
| a-20 | b+20 | c+20 |
| a+20 | b-20 | c-20 |
| a+20 | b-20 | c+20 |
| a+20 | b+20 | c-20 |
| a+20 | b+20 | c+20 |
Thus $3$ terms for each column make $E_x$ having $3 \times 8 = 24$ terms.
Part II : Finding divergence of first three terms of electric field
Let :
$\xi=a-20$
$\eta=b-20$
$\zeta=c-20$
As shown in the answer below by mathlove, $$\dfrac{\partial E_{x(1,2,3)}}{\partial a}=2\tan^{-1}\left( \dfrac{A}{C} \cdot \dfrac{E-D}{E-B} \right)$$
Therefore, for $\dfrac{\partial E_{x}}{\partial a}$, there is one term corresponding to each column in the above table. Therefore $\dfrac{\partial E_{x}}{\partial a}$ has $1 \times 8 = 8$ terms
Similarly, by repeating the steps from beginning in an analogous way, we find $\dfrac{\partial E_{y}}{\partial b}$ and $\dfrac{\partial E_{z}}{\partial c}$ each have 8 terms.
Thus $\vec{\nabla}.\vec{E}=\dfrac{\partial E_{x}}{\partial a}+\dfrac{\partial E_{y}}{\partial b}+\dfrac{\partial E_{z}}{\partial c}$ has $8+8+8 = 24$ terms.
How shall I proceed in getting the sum of entire $24$ terms as zero ?