Two spheres $a$ and $b$ both have point masses in their volumes. Each point mass is attracted to each point mass by an inverse-square law. The spheres both have a mass $M$ and $m$ and a moment of inertia $I$ and $i$.
The distance between the centres of $a$ and $b$ is $d$. The spherical coordinates $(r, \theta, \phi)$ of sphere $a$ are $(R, u, v)$ while those of sphere $b$ are $(r, w, x)$. Therefore the difference between $x$ components, $y$ components and $z$ components of $R$ and $r$ is $-R\sin(v)\cos(u)+ r\sin(x)\cos(w)+d$, $R\sin(v)\sin(u)-r\sin(w)(\sin x)$ and $R\cos v-r\cos x$ respectively. The distance from a point mass in $a$ to another in $b$ is: $$D={\sqrt{(-R\sin(v)\cos(u)+ r\sin(x)\cos(w)+d)^2+(R\sin(v)\sin(u)-r\sin(w)(\sin(x)^2+(R\cos(v)-r\cos(x))^2}}$$ The resultant force ($F$) is:- $$F=\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi}{K \over (-R\sin(v)\cos(u)+ r\sin(x)\cos(w)+d)^2+(R\sin(v)\sin(u)-r\sin(w)\sin(x))^2+(R\cos(v)-r\cos(x))^2}\,du\,dv\,dw\,dx$$
The component of this force in the $x$-direction is : $$f=\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi}{K(-R\sin(v)\cos(u)+ r\sin(x)\cos(w)+d)({\sqrt{(-R\sin(v)\cos(u)+ r\sin(x)\cos(w)+d)^2+(R\sin(v)\sin(u)-r\sin(w)(\sin(x)^2+(R\cos(v)-r\cos(x))^2}})\over (-R\sin(v)\cos(u)+ r\sin(x)\cos(w)+d)^2+(R\sin(v)\sin(u)-r\sin(w)\sin(x))^2+(R\cos(v)-r\cos(x))^2}\,du\,dv\,dw\,dx$$
where:- $$K= {GIi \over R^2r^2Mm}$$
If my working is correct, how can the integral above be solved to determine the value of f? Can the value of $f$ be obtained by any other suitable method?