I don't know, at which level you are studying Newtonian gravity, but it is very important in any case to have the right concepts about the different quantities in question.
First of all the gravitational force is a very fundamental concept, which cannot be derived somehow from simpler principles. It was Newton's ingenious insight to conclude from the empirical knowledge about celestial mechanics, i.e., the orbits of the planets around the Sun that there is a fundamental force acting between two point masses which is inversely proportional to the square of the distance of the bodies and their masses. It is always directed along the straight line connecting the two bodies and is always attractive. It is very important to keep in mind that a force is always a vector quantity.
For the gravitational force on a mass m_1 exerted by another mass m_2, located at the origin of your reference frame you obtain
\vec{F}_1=-\gamma \frac{m_1 m_2}{r^2} \frac{\vec{r}}{r},
where \vec{r} is the position of body 1.
Further the gravitational force from several bodies on body 1 is additive, i.e., the superposition principle holds. Thus for a point particle with mass m_1 the gravitational force due to a mass distribution \rho(\vec{x}) (where \rho is the mass density) you get
\vec{F}_1(\vec{r})=-\gamma m_1 \int_{\mathbb{R}^3} \mathrm{d}^3 \vec{x} \frac{\rho(\vec{x})(\vec{r}-\vec{x})}{|\vec{r}-\vec{x}|^3}.
The force is always proportional to the mass of the "test body" 1. Thus you can introduce the gravitational field a
\vec{G}(\vec{r})=\frac{\vec{F}_1}{m_1},
which is independent from the mass of the test body.
It is also immediately clear that the gravitational field is conservative, i.e., there exists a scalar potential for it such that the gravitational field is the negative gradient of this scalar potential:
\vec{G}(\vec{r})=-\vec{\nabla} V(\vec{r}).
For a single point mass in the origin of the coordinate system you have
\vec{G}(\vec{r})=-\gamma m_2 \frac{\vec{r}}{|\vec{r}|^3}.
It is easy to see that
V(\vec{r})=V(r)=-\frac{\gamma m_2}{r}, \quad r=|\vec{r}|.
That's easily seen by taking the derivatives carefully. Since the potential only depends on the distance we have
\vec{\nabla} V(r)=V'(r) \vec{\nabla} r.
Now
r=\sqrt{x^2+y^2+z^2} \; \Rightarrow \; \frac{\partial r}{\partial x}=\frac{x}{\sqrt{x^2+y^2+z^2}}=\frac{x}{r}
or analogously for the other two components
\vec{\nabla} r=\frac{\vec{r}}{r}.
Further we have
V'(r)=+\frac{\gamma m_2}{r^2}.
This finally leads to
\vec{G}=-\vec{\nabla} V(r)=-\frac{\gamma m_2}{r^2} \frac{\vec{r}}{r}=-\frac{\gamma m_2 \vec{r}}{r^3}.
With the more general case of a mass distribution you can use the superposition principle also for the potential. This leads to
V(\vec{r})=-\gamma \int_{\mathbb{R}^3} \mathrm{d}^3 \vec{x} \frac{\rho(\vec{x})}{|\vec{r}-\vec{x}|}.
You can as well bring this law to a local form by noting that this is the solution of the Poisson equation,
\Delta V(\vec{r})=4 \pi \gamma \rho(\vec{r}).
You can solve it by finding the Green's function of the Laplace operator.
This whole analysis ist prefectly analogous to the treatment of electrostatics. So you find most of the math about this subject well explained in textbooks on electromagnetism, e.g., Griffiths or Jackson.
