In the event you are interested in some methodology or a transparent way on how to solve this problem, here is an analytic approach.
Let the two spheres be $S_1$ and $S_2$ radius $R_1$ and $R_2$ at centres C1 and C2 respectively. The distance $C1C2 =d\ge R_1+R_2$. Work in spherical coordinates with the frame at C1:
1)Consider two point P1$( r_1\cos\phi_1\cos\theta_1, r_1\cos\phi_1\sin\theta_1, r_1\sin\phi_1)$
and P2$( d+ r_2\cos\phi_2\cos\theta_2, r_2\cos\phi_2\sin\theta_2, r_2\sin\phi_2)$ in the volume of $S1$ and $S2$ respectively.
2) Now write the square of the distance P1P2
$r_{12}^2=(d+ r_2\cos\phi_2\cos\theta_2- r_1\cos\phi_1\cos\theta_1)^2+(r_2\cos\phi_2\sin\theta_2-r_1\cos\phi_1\sin\theta_1)^2+ (r_1\sin\phi_1- r_2\sin\phi_2)^2$
3) Write the differential volume elements for the two spheres at the points P1 and P2 in spherical coordinates:
$dV_1=\sqrt {g_1}dr_1d\phi_1d\theta_1$ and $dV_2=\sqrt {g_2}dr_2d\phi_2d\theta_2$
where $g_1$ and $g_2$ are the Jacobian determinants for spherical coordinates (not difficult to find, or look it up in a vector analysis book, or use $dV=r^2\cos\phi drd\phi d\theta$). Hence the differential masses of these volume elements are
$dm_1=\rho_1(r_1)dV_1$ and $dm_2=\rho_2(r_2)dV_2$
Now you can put all these together and write the total magnitude of the force
$F=G\int_{V_1}\int_{V_2}\frac {\rho_1(r_1)\rho_2(r_2)dV_1dV_2}{r_{12}^2}$
and you need to do this integral with the transformations given above.
This is a very general formula giving the total force in the problem. One can see that if $d>>R_1+R_2$ this equation reduces to the gravitational attraction between two point masses at large distance $d$. One could end up with simpler integrals by using the symmetry in the line C1C2, and considering discs instead the general volume elements we have done in the above method