About negative energies: they set no problem:
On this context, only energy differences have significance. Negative energy appears because when you've made the integration, you've set one point where you set your energy to 0. In this case, you have chosen that $PE_1 = 0$ for $r = \infty$. If you've set $PE_1 = 1000$ at $r = \infty$, the energy was positive for some r.
However, the minus sign is important, as it is telling you that the test particle is losing potential energy when moving to $r = 0$, this is true because it is accelerating, causing an increase in $KE$:
let's calculate the $\Delta PE_1$ for a particle moving in direction of $r = 0$: $r_i = 10$ and $r_f = 1$:
$\Delta PE_1 = PE_f - PE_i = Gm(-1 - (-0.1)) = -Gm\times0.9 < 0$
as expected: we lose $PE$ and win $KE$.
Second bullet: yes, you are right. However, it is only true IF they are point particles: has they normally have a definite radius, they collide when $r = r_1 + r_2$, causing an elastic or inelastic collision.
Third bullet: you are right with $PE_2 = mgh$, however, again, you are choosing a given referential: you are assuming $PE_2 = 0$ for $y = 0$, which, on the previous notation, means that you were setting $PE_1 = 0$ for $ r = r_{earth}$.
The most important difference now is that you are saying that an increase in h is moving farther in r (if you are higher, you are farther from the Earth center).
By making the analogy to the previous problem, imagine you want to obtain the $\Delta PE_2$. In this case, you begin at $h_i = 10$ and you want move to $h_f = 1$ (moving in direction to Earth center, like $\Delta PE_1$:
$\Delta PE_2 = PE_{f} - PE_{i} = 1mg - 10mg = -9mg < 0$.
As expected, because we are falling, we are losing $PE$ and winning $KE$, the same result has $PE_1$
Fourth bullet: they both represent the same thing. The difference is that $gh$ is the first term in the Taylor series of the expansion of $PE_1$ near $r = r_{Earth}$. As exercise, try to expand $PE_1(r)$ in a taylor series, and show that the linear term is:
$PE_1 = a + \frac{Gm(r-r_{earth})}{r_{earth}^2}$.
Them numerically calculate $Gm/r_{earth}^2$ (remember that $m=m_{earth}$). If you haven't made this already, I guess you will be surprised.
So, from what I understood, your logic is totally correct, apart from two key points:
energy is defined apart of a constant value.
in the $PE_1$, increase r means decrease $1/r$, which means increase $PE_2 = -Gm/r$. In $PE_2$, increase h means increase $PE_2=mgh$.
The short answer is "observations".
In the case of the gravitational law, the orbits of the planets around the sun, the moon around the earth fit mathematically a force with an inverse square law for the distance. An inverse law does not.
In the case of electricity this article points out the observational history:
Early investigators who suspected that the electrical force diminished with distance as the gravitational force did (i.e., as the inverse square of the distance) included Daniel Bernoulli 1 and Alessandro Volta, both of whom measured the force between plates of a capacitor, and Aepinus who supposed the inverse-square law in 1758.
and then others took it from there to end up with the comprehensive publications of Coulomb, based on measurements.
Finally, in 1785, the French physicist Charles Augustin de Coulomb published his first three reports of electricity and magnetism where he stated his law. This publication was essential to the development of the theory of electromagnetism.[8] He used a torsion balance to study the repulsion and attraction forces of charged particles and determined that the magnitude of the electric force between two point charges is directly proportional to the product of the charges and inversely proportional to the square of the distance between them.
Best Answer
The formula:
$$ \Delta U = mgh $$
is an approximation that applies when the distance $h$ is small enough that changes in $g$ can be ignored. As you say, the expression for $U$ is:
$$ U= -\frac{G M m}{r} $$
So the change when moving a distance $h$ upwards is:
$$ \Delta U = \frac{GMm}{r} - \frac{GMm}{r + h} $$
We rearrange this to get:
$$\begin{align} \Delta U &= GMm \left( \frac{1}{r} - \frac{1}{r + h} \right) \\ &= GMm \frac{h}{r^2 + rh} \\ &= \frac{GM}{r^2} m \frac{h}{1 + h/r} \\ &\approx \frac{GM}{r^2} m h \end{align}$$
where the last approximation is because $h \ll r$ so $1 + h/r \approx 1$. And since $GM/r^2$ is just the gravitational acceleration $g$ at a distance $r$, we get:
$$ \Delta U = g m h $$