The quadratic potential is the simplest possible model for a bond. You can derive it by considering the Taylor expansion of the potential around the natural bond length
$V(r - r_0) = V(r_0) + \frac{d V(r_0)}{d r} (r - r_0) + \frac{1}{2} \frac{d^2 V}{dr^2} (r - r_0)^2$
The constant term can be set to 0 since it does not contribute to the force and just sets the 0 point of your energy scale. The linear term is 0, because you are at a stationary point [1].
That leaves you with
$V(r - r_0) = \frac{1}{2} \frac{d^2 V}{dr^2} (r - r_0)^2$
If we require the second derivative to be constant, we have recovered the harmonic potential.
As @Charlie Crown illustrated, this force resulting from this potential does not go to 0 at infinity, while the Morse potential does. You can of course take a polynomial of higher order than two, but not every order is suitable. A third order polynomial results in a potential that (typically) goes to negative infinity at large $r$, so instead a quartic potential is some times used. It has the advantage of being slightly "wider" than the quadratic one. That said, a completely unrelated reason why none of these can simulate bond breaking/formation is that the implementation requires explicit declaration of which atoms should interact via the stretching potential.
Still, at large $r$ both differ significantly from the Morse potential. Why then is the Morse potential not used?
The restoring force for large $r$ is very low in case of the Morse potential, hence it takes longer for the bond length to return to the equilibrium position. The quadratic potential describes the potential well for displacements close to equilibrium and for moderate temperatures, this is the part of the potential you care about.
Obviously that still leaves the question of how to simulate bond breaking in a force field. ReaxFF assumes that the bond order of a pair of atoms can be determined from the interatomic distance alone.

(qualitative recreation from [2])
The sigma, pi and double pi bonds contribute increasingly to the overall bond order (max individual bond order is 1) as the atoms get closer together. For simplicity I am leaving out the corrections made to the overall bond order necessitated by overcoordination.
The bond stretching potential takes the form of a modified Morse potential
$E_{Bond} = -D_e \cdot BO_{ij} \cdot \exp(p \cdot (1 - BO_{ij}^p))$
where $p$ is a bond specific parameter[2].
References:
[1]: Frank Jensen, Introduction to Computational Chemistry Chap. 2
[2]: J. Phys. Chem. A 2001, 105, 9396-9409