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.