Some people do:
In this paper there is a system coupled to a bath of Morse oscillators rather than a bath of harmonic oscillators, but it is not exactly solvable, they used a numerical approach called mctdh. When it is said that the Morse potential is "exactly solvable", what it means is that you can solve the vibrational Schroedinger equation for a Morse potential analytically (the end result is a formula for the vibrational energies and wavefunctions of the system). What about if there is a spin coupled to 500 Morse oscillators, do we have an exact solution for the overall energies and wavefunctions of this very complicated system? We do not even have an exact solution in the case where the oscillators are harmonic oscillators; this would be called the "spin-boson problem" and it's not exactly solvable except in specific cases. Finally, even if there was only one Morse oscillator, the solutions that you see in that Wikipedia page you linked, are not very simple or easy to use: For example the exact solutions for the quantum harmonic oscillator do not involve generalized Laguerre polynomials but for the Morse oscillator they do.
Most people don't:
In solid state physics, many people model a system coupled to phonons using models like the spin-boson model or generalizations of it. For example in this paper of mine we study how a qubit would undergo decoherence if a quantum computer were to be made with GaAs quantum dots. The qubit can be defined as follows: 0 = absence of an exciton, 1 = presence of an exciton, so it's a 2-level system, but it's coupled to all of the vibrations of the semiconductor lattice in which it lies. 2-level systems can be considered as "spin-1/2 particles", so what we have is a spin interacting with a bunch of vibrations. These vibrations are approximated to be harmonic oscillators, so we simply have the aforementioned "spin-boson problem" which has been studied for several decades (and still cannot be solved exactly most of the time). Now there's several reasons why we chose this simple spin-boson model with harmonic oscillators rather than using a Morse potential:
- Neither case (harmonic oscillators or Morse oscillators) is exactly solvable when you have not only the oscillators (nuclear vibrations) but also the spin (or electronic/excitonic degrees of freedom), but at least for the harmonic case we have simple analytic expressions for things such as the Feynman-Vernon influence functional which describes the influence of the vibrations on the spin; for Morse oscillators we do not have such a simple influence functional. So the calculations are much easier in the harmonic case compared to the Morse oscillator case.
- A Morse oscillator is actually harmonic at the very bottom, and only starts to deviate from a harmonic oscillator for much larger internuclear distances as the system starts to dissociate. If we were to model a qubit in a GaAs semiconductor at such a high temperature that the semiconductor is in the midst of breaking apart (the Ga and As atom are dissociating from one another), we may wish to use an anharmonic potential, but this is rarely (or never!) the case. Consider this: Which vibrational levels of these oscillators actually have a significant population in your system? Surely not the ones at the top of this Morse potential (see v=8 and 9 and observe that their outer turning points are at internuclear distances of around r=12 and 14 which is about triple the internuclear distance at equilibrium; do you picture your solid-state lattice having its internuclear distances so far away from equilibrium in the normal scenario in your interest?).
- For this particular system, a spectral distribution function (which tells you how strongly each oscillator couples to the "spin") was determined from fitting the dynamics of the spin-boson model, and the empirically obtained parameters of the spectral function matched almost exactly the values of those parameters obtained from first-principles calculations, and the chances of this being pure luck are so extremely low that we believe that the simple spin-boson model with harmonic oscillators is an excellent approximation.
- If we were not in a case where we know the spin-boson model with harmonic oscillators is an excellent approximation (as described in the last bullet point), we can stop and think about all the other multitude of approximations we are making (maybe the ignoring of spin-orbit coupling, maybe the ignoring of the difference between relativistic mass and non-relativistic mass, maybe the use of the Born-Oppenheimer separation between the electronic/excitonic and nuclear/phononic/bosonic degrees of freedom, or maybe we don't have the exact parameters for the exact solid in question so we are just assuming we can use the parameters from a similar solid that has been studied in more detail, etc.), and then realize that there's just so many approximations going on that Jon von Neumann's quote applies: "Why be precise when we don't know what we're talking about?" Will you use quadruple precision (33-36 digits) for solving a differential equation telling you what weather it will be tomorrow, when the coefficients in the differential equation have enormous error bars? Then don't use an anharmonic bath if you don't have a fairly precise description of all other relevant information (relativistic effects, spin-orbit coupling, etc.).
Conclusion: The bottom line is that you would be making your life harder (more difficult equations that can only be solved less efficiently or with less accuracy with the same resources), and most of the time you would not be gaining a better understanding of the relevant physics.
