Maximum Gravity Around a Unit Sphere Most simulations that involve planetary gravity use Newton's law of universal gravitation and treat planets like point-masses. This is very accurate at large distances, and fairly accurate all the time, but introduces small errors when close to the surface. I wanted to find the distance from a planetary body at which gravity is at a maximum. My "planet" is a homogeneous unit sphere, and I know that when r=0, gravity is zero, and at large r, gravity is roughly 1/r^2. At some distance near the surface, perhaps on the surface, there is a maximum. This is the equation that I came up with to find this maximum:

where r is the distance from the center of the body along the x axis. I think this will be accurate as long as r>=1, but Mathematica can't evaluate the integral. My question: is this equation correct, and how can I solve it.
 A: You assume the force to be
$$ F = \frac{1}{(x-r)^2 + y^2 + z^2} $$
but you have to take into account that different parts of the sphere will pull you in different directions. So you have to multiply with the normalized vector:
$$\begin{pmatrix} r-x \\ -y \\ -z \end{pmatrix}\frac{1}{\sqrt{(x-r)^2 + y^2 + z^2}}$$
giving in $x$ direction:
$$ F = \frac{r-x}{\left((x-r)^2 + y^2 + z^2\right)^{3/2}} $$
I don't think you can get Mathematica to solve this since the result is pretty singular. You could just imploy the shell theorem which tells you that if you are outside the sphere the gravity is exactly as for a point mass. If you are inside the sphere it is the same, except that you need to consider just the mass of the sphere that is inside the radius of $r$.
Outside the sphere $(r>1)$ the total force is:
$$ F = \frac{M_0}{r^2} $$
The mass inside a radius $r$ is proportional to $r^3$:
$$ M(r) = M_0 r^3 $$
Inside the sphere $(r<1)$ we get:
$$ F = \frac{M(r)}{r^2} = M_0 r $$
For $r=1$ the two agree giving you the maximum there.
A: I found this question a little confusing at first, until I realized with the help of the comments that the body as described in the problem is apparently a sphere of radius $1$ and center at the point $(x, y, z) = (r, 0, 0)$, having constant unit density, and the integrals given express the field exterior to the body.  At least I think that's correct.  If it is, it may help explain why Mathematica choked on the integrals given in the question:  there is no position $(r, 0, 0)$ at which the integrals are stationary with respect to $r$, since magnitude of the force decreases monotonically to 0 as the distance to the body grows, and its derivative (with respect to the distance to the body's center) never vanishes:  there is no stationary solution.  But the following argument shows how the problem may be solved.  I guess the procedure presented below is related to the so-called shell theorem, but I did not know it by that name, and worked out the derivations as I went along from first principles, aided no doubt by memory.  But I certainly lay no claim to inventing this argument. 
The easiest way I know to solve this problem is to exploit spherical symmetry, Gauss's law, and the fact that the gravitational potential $\phi$ satisfies a Poisson equation
$\nabla^2 \phi = 4 \pi G \rho, \tag{1}$
where $\rho$ is the mass per unit volume and $G$ is Newton's universal constant of gravitation; see this Widipedia entry  Then the gravitational force $\mathbf F$ experienced by a test particle, is, per unit mass, given by
$\mathbf F = -\nabla \phi. \tag{2}$
If we envision a spherical body of radius $R$ and uniform mass density $\rho_0$ centered at $r = 0$ we can use spherical symmetry to evaluate the surface integral of the force $\mathbf F$  on the spherical surface $\partial S_r$ of the sphere $S_r$ of radius $r$, again centered at $r = 0$.  We have 
$\int_{\partial S_r} \mathbf F \cdot \mathbf n dA = \mathbf F \cdot \mathbf n \int _{\partial S_r} dA, \tag{3}$
since spherical symmetry implies that $\mathbf F \cdot \mathbf n$ in fact takes a constant value on $\partial S_r$.  Here we take $\mathbf n$ to be the outward pointing unit normal vector on $\partial S_r$; furthermore we use the fact that $\mathbf F$ must everywhere be collinear with $\mathbf n$.  This follows from (2) and the assumption that $\phi$ is itself spherically symmetric; as such, it is constant on $\partial S_r$ for any $r > 0$; thus, any component of $\nabla \phi$ in a direction tangent to $\partial S_r$ must vanish, leaving $\mathbf F = -\nabla \phi$ with only a radial component; that is, one collinear with $\mathbf n$.  We can thus write $\mathbf F = \mathbf F_r \mathbf n$, where $\mathbf F_r$ is constant on $\partial S_r$, again by spherical symmetry.  Now since $\mathbf F \cdot \mathbf n = \mathbf F_r \mathbf n \cdot \mathbf n = \mathbf F_r$, we can bring it out of the integral on the left-hand side of (3), and that is in fact how (3) is proved.  Now the right-hand side of (3) may be expressed thusly, introducing into it the notation $\mathbf F_r$ for $\mathbf F \cdot \mathbf n$:
$\mathbf F_r \int_{\partial S_r} dA = 4 \pi r^2 \mathbf F_r, \tag{4}$
since $\int_{\partial S_r} dA$ is simply the area of the sphere of radius $r$.  Inserting (4) into (3) and re-arranging, we arrive at a formula for $\mathbf F_r$:
$\mathbf F_r = \frac{1}{4 \pi r^2}\int_{\partial S_r} \mathbf F \cdot \mathbf n dA, \tag{5}$
and use Gauss's theorem to evaluate $\int_{\partial S_r} \mathbf F \cdot \mathbf n dA$ in the two cases $r \le R$ and $r > R$; the theorem tells us that
$\int_{\partial S_r} \mathbf F \cdot \mathbf n dA = \int_{S_r} \nabla \cdot \mathbf F dV, \tag{6}$
where $dV$ is the volume element, and from (1) and (2) we see that
$\nabla \cdot \mathbf F = - \nabla^2 \phi = -4 \pi G \rho. \tag{7}$
Combining (5), (6), and (7) yields
$\mathbf F_r = -\frac{G}{r^2} \int_{S_r} \rho dV. \tag{8}$
We evaluate the integral occurring in (8) in each of the two cases at hand.  For $r \le R$, $\rho = \rho_0$, whence
$\mathbf F_r = -\frac{G}{r^2} \int_{S_r} \rho_0 dV = -\frac{G \rho_0}{r^2} (\frac{4}{3} \pi r^3) =  -\frac{4}{3} G \rho_0 \pi r; \tag{9}$
if, on the other hand, $r > R$, we obtain
$\int_{S_r} \rho dV = \frac{4}{3}\pi \rho_0 R^3, \tag{10}$
since the the net mass inside a sphere of radius $r$ is in this case in fact all contained in the sphere of radius $R$, which is "full", as it were; the rest of the space is empty for $r > R$.  $\frac{4}{3}\pi \rho_0 R^3$ is in fact the net mass of the sphere $S_R$; it is seen that $\mathbf F_r$ is continuous in $r$:  the inside ande outside formulas agree at $r = R$; furthermore, $\mathbf F$ is always directed inward, in the direction of $-\mathbf n$.  Inserting (10) into (8) we arrive at
$\mathbf F_r = -\frac{4}{3} G \rho_0 \pi \frac{R^3}{r^2}. \tag{11}$
Combining (9) and (11) we can say:
The gravitational force of a spherical body of uniform mass density centered at $r = 0$ takes the form
$\mathbf F = -\frac{4}{3}G \rho_0 \pi f_R(r) \mathbf n, \tag{12}$
where $f_R(r)$ is the function from $[0, \Bbb R) \to \Bbb R$ defined by
$f_R(r) = r, \, 0 \le r \le R, \tag{13}$
$f_R(r) = \frac{R^3}{r^2}, \, R \le r. \tag{14}$
As such we can say that the magnitude $\mathbf F_r$ of the force $\mathbf F$,
$\vert \mathbf F_r \vert = \frac{4}{3}G \rho_0 \pi f_R(r), \tag{15}$
increases from zero linearly with $r$ from $r = 0$ until $r = R$, at which point it decreases according to $r^{-2}$.  The maximum field strength is reached at $r = R$, with a value
$\vert \mathbf F_r \vert = \frac{4}{3} G \rho_0 \pi R. \tag{16}$
Hope this helps.  Cheers, 
and as always,
Fiat Lux!!!
