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.
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.