Can the Schrodinger equation describe planetary motion? I was asked on an exam whether the Schrodinger equation can be used to describe planetary motion and my answer was "No, because the solutions are wavefunctions which give probabilities but everything can be exactly measured for large objects."
Then I read this article which suggested it should be possible. What would the Hamiltonian be, and how do we get definite results instead of probabilities?
 A: Yes, you can construct classical orbits from the Schrodinger equation, as long as you take the right limit. For example, consider the hydrogen atom. While the lower energy levels, like the $1s$ or $2p$ look nothing like classical orbits, you may construct wavefunctions that do by superposing solutions with high $n$ appropriately. Such solutions obey the Schrodinger equation with Hamiltonian $H = p^2/2m - e^2/r$, but have a sharp peak, which orbits the nucleus in a circular or elliptical trajectory. For very high $n$ the peak can get extremely sharp, so the position is definite for all practical purposes.
Heuristically this works because, given a state that is a superposition of states with different $\ell$'s and $m$'s but the same $n$, the coefficients of the $|n, \ell, m \rangle$ states inside are essentially a discrete Fourier transform of the position space wavefunction, with $O(n^2)$ entries. For higher $n$, you can use this to build sharper and sharper peaks. This reasoning is not even necessary; we know it has to work because everything is quantum, so there must be a way to reproduce classical results within quantum mechanics.
By treating a planet as a single particle, the same reasoning holds. However, as pointed out by other answers, this isn't the full story, because a planet is much more complicated than an electron. In fact this complication is essential, because if a planet and star were single particles in perfectly empty space, there's no particular reason the planet would end up in one of these classical-looking states with a sharply peaked position. For a real planet this is a consequence of decoherence: superpositions that don't have sharply peaked positions are not stable against interaction with the environment. That's how the classical world emerges.
I should also note that the link you gave is not an example of this. That article is about an astrophysical quantity that happens to be described by an equation with the same form as the Schrodinger equation, but absolutely nothing quantum is going on. It's just a coincidence, which arises because there aren't that many different low-order PDEs. If there is any reason at all, it's simply that we look for simple equations on both microscopic and macroscopic scales. 
A: Yes, no and no.
Yes in the sense that insofar everything is described by quantum mechanics it will evolve according to the Schrödinger equation.
No in the sense that planets are composed by a lot of particles, and the equation will have $3N$ dimensions describing the joint wave-function. It will also have a lot of complicated interaction terms (and since we are talking about fermions it will actually be the Dirac equation). So it will not be an analogue of the simple electron around proton spherically symmetric equation textbooks solve for hydrogen, but a huge monster version that is likely far from what the exam considered. We know empirically that planetary motion nicely averages together into the classical orbit, but why that happens is separate from the Schrödinger equation framework. 
Finally, no in the sense that proper orbits require general relativity and hence are incompatible with quantum mechanics. One can  try to solve the Schrödinger equation in a Schwartzschild spacetime (i.e. assume the planet has no effect on the central body) as a semiclassical model: this produces a somewhat hydrogen-atom like result for "particle" planets with no components. But this is obviously only an approximation since there will be a backreaction on the central body. 
A: 
I was asked on an exam whether the Schrodinger equation can be used to describe planetary motion and my answer was "No, because the solutions are wavefunctions which give probabilities but everything can be exactly measured for large objects."

There are some problems with your answer. First, wavefunctions do not just describe probabilities. In some situations the square modulus of the wave function obeys the rules of probability, but in many experiments it does not:
https://arxiv.org/abs/math/9911150
Second, it's not clear what you mean by "everything can be measured exactly for large objects". Quantum mechanics sez that uncertainties in position and momentum obey the following inequality $\Delta x\Delta p \geq \hbar$. I don't think astronomers measure planets with a high enough accuracy to measure that small an effect.
Third, quantum mechanics has been applied to planetary dynamics:
https://arxiv.org/abs/quant-ph/9612037
Quantum mechanics predicts that large objects like planets won't undergo interference due to decoherence caused by interacting with other systems, such as light incident on the planet in question.

Then I read this article which suggested it should be possible. What would the Hamiltonian be, and how do we get definite results instead of probabilities?

The article you linked to is not relevant to the issue you raised since it describes the Schrodinger equation as an approximate description for a load of particles interacting, not as a description of a single particle or object as it is used in quantum mechanics.
A: It's very common in physics that the same differential equation describes different systems. For example consider the equation:
$$ \frac{dy}{dt} = -k\,y $$
This describes decay of radioactive substances, the rate of water flow through a leak at the bottom of a barrel, and probably the number of live brain cells in my brain. But even though it's the same equation in these examples the physical meaning of $y$ and the constant $k$ are very different.
Now, the work you mention is described in Schrödinger Evolution of Self-Gravitating Disks by Konstantin Batygin and the details are described in that paper. Batygin constructs a function $\psi$ that is a rather abstract function constructed from variables that describe the motion, and he finds that $\psi$ obeys the Schrodinger equation. But it is physically completely different. His $\psi$ is not a wavefunction and is unrelated to probabilities.
A: The Schrödinger equation occurring in the article is a classical equation of motion, that effectively describes some continuum mechanics problem. In this sense, the disks are not governed by the Schrödinger equation with its quantum interpretation. Similarly, Schrödinger-like equations occurs as classical wave equations in may areas (e.g. water waves).
However, the actual Schrödinger equation can and does describe the orbital motion of planets. Since Newtonian mechanics must be a limit of quantum mechanics, if we want to accept quantum mechanics as a model describing our world (after all, the new theory has to explain why the old one was so successful).
The simplest way to see, that this is the case is to study the Ehrenfest theorem, which tells us how expectation values evolve. The derivation from the Schrödinger equation and its conjugate is quite straight forward:
\begin{align*}
 i\hbar \partial_t \left|\psi\right> &= H\left|\psi\right> &
 -i\hbar \partial_t \left<\psi\right| &= \left<\psi\right| H
\end{align*}
\begin{align*}
\partial_t \left<\psi \middle| A \middle| \psi \right> &= (\partial_t \left<\psi\right|) A \left|\psi\right> + \left<\psi\right| (\partial_t A) \left|\psi\right> + \left<\psi\right|A \partial_t \left|\psi\right> \\
&= \frac i \hbar \left<\psi\right| HA \left|\psi\right> + \left< \psi \right| (\partial_t A) \left|\psi\right> - \frac i \hbar \left<\psi\right|AH\left|\psi\right>
\end{align*}
$$ \partial_t \left< A \right> = \frac i \hbar \left<\left[H, A\right]\right> + \left<\partial_t A\right> .$$
If we now take a $n$-particle Hamiltonian:
$$ H = \sum_i \frac{p_i^2}{2m_i} + V(\vec r_1, \ldots, \vec r_n)$$
and write down the equations for the expectation values for $\vec r_i$ and $\vec p_i$ we get the equations:
\begin{align*}
  \partial_t \left<\vec p_i \right> &= \frac i \hbar \left<[H, \vec p_i]\right> = \left< \nabla_i V(\vec r_1, \ldots,  \vec r_n) \right> \\
  \partial_t \left<\vec r_i \right> &= \frac i \hbar \left<[H, \vec r_i]\right> = \frac i \hbar \left<\left[\frac{p_i^2}{2m_i}, \vec r_i\right]\right> = \frac {\left<\vec p_i\right>} {m_i}. \\
\end{align*}
Those are almost the classical equations of motion for $\left<\vec r_i\right>$ and $\left<\vec p_i\right>$. The only difference is, that the force $\nabla V(\vec r_1, \ldots, \vec r_n)$ is not evaluated at the average position, but averaged over state. This, however, does not matter if our state is a very sharp wave packet compared to the length scale on which $\nabla_i V$ varies.
This shows how quantum mechanics describes orbital motion, since it describes classical mechanics in the limit of sharply concentrated wave packets in the ($\vec p$, $\vec r$) plane, so measuring will have a very small quantum mechanical uncertainty. The limits put on the sharpness of the peak by the uncertainty relation are negligible for a planet${}^1$ which has immense mass, so that even for small velocities $\vec p$ gets very large compared to $\hbar$, so that the position uncertainty can also be very small.
If we have some bound system made up of $n$ particles, we can do the same thing as in classical mechanics and transform the coordinates to get an equation for the centre of mass and equations for the relative motion of the constituents. That, in turn, means that we can write down the Schrödinger equation for the centre of mass of a planet, just like we can write down equations of motion for the centre of mass in classical mechanics. (Note: that 
potential that acts on the centre of mass is exactly the same as if the planets were point masses requires spherical symmetry of the orbiting bodies and is related to Newtons shell theorem – but it is true for more general bodies as the fist order of a multipole expansion).
Using the idea of the last paragraph, we can write down the Hamiltonian for an orbiting planet. It has the same form as the Hamiltonian for a hydrogen atom:
$$ H = \frac{p_1^2}{2m_1} + \frac{p_2^2}{2m_2} + G\frac{m_1 m_2}{\left|\vec r_1 - \vec r_2\right|}, $$
Of course it gets more complicated if we include several masses, but the presented Hamiltonian has the convenient property, that the solutions are known, so we can study exact solutions.
Now, we can even go further and pose the question how the classical orbits arise as time dependent bound state solution$^2$ of the Schrödinger equation (which are, after all, proportional to spherical harmonics, so smeared over the whole orbit). To do this we have to construct the sharp wave packets (the discussion follows the one at the end (p. 136-138) of Chapter 6.3 of Franz Schwabl: Quantum Mechanics. Fourth Edition, Springer (2007)). We restrict the discussion to circular orbits, construction of sharply peaked solutions on elliptical orbits is much more difficult. In analogy to the way Gaussian wave packets are derived for free particles, we superimpose eigenstates with large quantum numbers $n$. Further, we choose only those components with maximal angular momentum $l = n-1$ (since the classical circular orbits have the largest angular momentum for a given energy) and maximal magnetic quantum number $m = l$ (since those are the states that are most localized around the central plane of rotation). This gives a form
\begin{align*}
 \psi(r, \vartheta, \varphi, t) &= \sum_n c_n \psi_{n,n-1,n-1}(r, \vartheta, \varphi) e^{-iE_nt/\hbar}
\end{align*}
The $c_n$ are chosen to be centred around some large $n_0$ and decay on a width small compared to $n_0$. We can now write $n = n_0 + \varepsilon$ and develop in the small parameter $\varepsilon/n_0$
$$ \psi(r, \vartheta, \varphi, t)  = \sum_n c_n \frac{1}{\sqrt \pi n! n^n a^{3/2}} \left(-\frac r a \sin(\vartheta)\right)^{n-1} e^{-r/na} e^{i(n_0+\varepsilon)\varphi + i t \frac{\left|E_0\right|}{\hbar} \left(1/n_0^2 - 2\varepsilon/n_0^3\right)} $$
This is obviously has a sharp peak around $\theta = \pi/2$ for large $n_0$. In the $\theta = \pi/2$ plane, we can write (the normalization and constant phase is absorbed in the new development coefficients $c_n'$):
$$ \psi(r, \vartheta=\pi/2, \varphi, t) \propto \sum_n c_n' r^{n-1} e^{-r/na} e^{i(n-n_0) \big(\varphi - (2\left|E_0\right|/\hbar n_0^3)t\big)}$$
If the $c_n'$ are chosen appropriately, this will show a time dependence $\psi(r, \pi/2, \varphi, t) = f(r)g(\varphi-\omega t)$ with $\omega = 2\left|E_0\right|/\hbar n_0^3$. The radial distribution is similarly sharply peaked (it is a bit more work to show this, the easiest is to consider that relative difference of the position expectation values $\frac{\left<r\right>_{n,n-1,n-1} - \left<r\right>_{n+1,n,n}}{\left<r\right>_{n,n-1,n-1}} \approx 2/n \to 0$ for large $n$, but it is a bit more work to show that the radial uncertainty also gets small as an absolute number).
The nice thing about these wave packets is that, unlike Gaussian wave packets in free space, their uncertainties do not grow unboundedly. So we do not even have to worry that the position of our planet will be smeared all over the solar system after a few million years. (This of course, can not happen for macroscopic objects, but discussing this would lead too far away from the question.)
Another remark is, that such quasi-classical orbits can also be prepared for atoms. Such highly excited atoms are called Rydberg atoms.

${}^1$ If we want to realize an uncertainty of $\Delta x = 10^{-30}\,\mathrm{m}$ of the earth's centre the velocity's uncertainty will be bounded by $\Delta v \ge \frac{\hbar}{2 m_E \Delta x} = 8.5 \cdot 10^{-30}\,\mathrm{\frac m s}$.
${}^2$ The bound state is important here in so far, as that it is simple to use a Gaussian wave packet and the Ehrenfest apparatus derived above. The problem is, that (an exponentially small) part of the probability will escape to infinity then, because the probability for arbitrarily high momentum is not zero. So the solution will not be superposition of bound state solutions but include some coefficients in the continuous spectrum.
