Here is my attempt at derivation:
First you must find a vector function for the gravitational force.
By the inverse square law, the magnitude of gravitational force between two bodies of mass $m$ and $M$ of distance $r$ apart is: $$G \frac{M m}{r^2}$$
The direction of this force points towards the other body. If you let $\vec{r}$ be the position vector from the other body towards you, then $\frac{-\vec{r}}{\|\vec{r}\|}$ gives you a unit radial vector pointing in the direction of the other body. This can be scaled by the magnitude to give the force vector as:
$$F(\vec{r}) = G \frac{M m}{\|\vec{r}\|^2} * \frac{-\vec{r}}{\|\vec{r}\|} = -G \frac{M m}{\|\vec{r}\|^3}\vec{r}$$
Now, potential energy is defined as
$$ U = -W = -\int_C F(\vec{r}) \cdot d\vec{r}$$
The path along which we take the line integral needs to be from an infinite distance away (which we set as our reference point with zero potential) to our current position.
The radial displacement vector $\vec{r}$ can be broken up into a component-wise displacement vector:
$$\vec{r} = x\hat{i} + y\hat{j} + z\hat{k}$$
The differential $dr$ is then
$$d\vec{r} = dx\hat{i} + dy\hat{j} + dz\hat{k}$$
Substituting into the expression for the dot product, we get that:
$$ F(\vec{r}) \cdot d\vec{r} = -G \frac{M m}{\left(\sqrt{x^2 + y^2 + z^2}\right)^3}(x\hat{i} + y\hat{j} + z\hat{k}) \cdot (dx\hat{i} + dy\hat{j} + dz\hat{k})$$
This gives
$$-\frac{G m M x}{\left(x^2+y^2+z^2\right)^{3/2}}dx -\frac{G m M y}{\left(x^2+y^2+z^2\right)^{3/2}}dy-\frac{G m M z}{\left(x^2+y^2+z^2\right)^{3/2}}dz$$
This can then be integrated term by term. The limits of integration of each should vary from $\infty$ to the current position ($x$, $y$, or $z$ since gravity is a conservative force (which can be verified mathematically by checking if the y-partial of the x component and the x-partial of the y component are equal) and thus path you take when computing work done does not matter -- it depends only on the initial and final positions.
This gives
$$\int_{\infty }^x -\frac{G m M x}{\left(x^2+y^2+z^2\right)^{3/2}} \, dx + \int_{\infty }^y -\frac{G m M y}{\left(x^2+y^2+z^2\right)^{3/2}} \, dy + \int_{\infty }^z -\frac{G m M z}{\left(x^2+y^2+z^2\right)^{3/2}} \, dz$$
This evaluates to
$$\frac{G m M}{\sqrt{x^2+y^2+z^2}} + \frac{G m M}{\sqrt{x^2+y^2+z^2}} + \frac{G m M}{\sqrt{x^2+y^2+z^2}} = \frac{3 G m M}{\|r\|}$$$$\frac{G m M}{\sqrt{x^2+y^2+z^2}} + \frac{G m M}{\sqrt{x^2+y^2+z^2}} + \frac{G m M}{\sqrt{x^2+y^2+z^2}} = \frac{3 G m M}{\|\vec{r}\|}$$
Then taking potential as the negative of the work done we get
$$U = -\frac{3 G m M}{\|r\|}$$$$U = -\frac{3 G m M}{\|\vec{r}\|}$$ However, this is clearly incorrect as there is a factor of three that should not be there.
Did I integrate incorrectly or am I missing something else?
