Does classical electromagnetism really predict the instability of atoms? I will try to give a concise summary of what I wrote below. I understand that it is very long and apologize if I am wasting your time.
I used the Liénard-Wiechert potential and the Lorentz force formula to derive equations of motion for charged particles that do not involve the $\mathbf{E}$ and $\mathbf{B}$ fields. I use this in a two-body problem setup with two charges, and want to see whether the resulting equations result in spiraling orbits.
More specifically, I plug in the formulas for $\mathbf{E}_2$ and $\mathbf{B}_2$ given by Liénard-Wiechert into the Lorentz force formula $\mathbf{F}_{12}(\mathbf{r}_1(t), t) = q_1(\mathbf{E}_2(\mathbf{r}_1(t), t) + \mathbf{\dot{r}}_1(t)\times\mathbf{B}_2(\mathbf{r}_1(t), t))$ and use the classical equation of motion $m_1\mathbf{\ddot{r}}_1(t) = \mathbf{F}_{12}(\mathbf{r}_1(t), t)$. A second equation is obtained for particle 2 using the expressions for $\mathbf{E}_1$ and $\mathbf{B}_1$ instead.
The full equations are complicated so I take a limit as the mass of the nucleus goes to infinity, and (assuming my calculations are correct), I obtain a simpler equation of motion formally identical to that of planetary motion around the sun. This predicts elliptical trajectories for most initial conditions.
The concern is that the resulting motion is stable, and this seems to contradict what the usual accounts say about classical electromagnetism failing to explain atomic orbit stability.
If you want more detail, you can also read what I wrote below.

Suppose for a moment that we work classically and use the planetary model for the atom (hydrogen, for simplicity), as a positively charged nucleus with a negatively charged electron orbiting it in a similar way that planets orbit around the Sun.
The usual story goes that due to the electron orbiting around the nucleus, it undergoes acceleration which will cause the electron to radiate electromagnetic waves. The energy of this radiation is considered to be taken away over time from the electron's total energy via the Larmor formula, and therefore the model predicts a collapse of the electron's orbit over a short period of time because the radius of the orbit must decrease to compensate for the diminished energy of the electron.
At the risk of sounding ridiculous to more knowledgeable people, I would like to challenge this assumption with the following considerations (only as a way to clarify my own understanding of the problem). It appears to me that this problem is arising only because we consider the electromagnetic field to have an existence independent of the charges generating it, and carrying energy and momentum on its own.
But it seems possible to describe the fundamentals of (classical) electromagnetism without resorting to the concept of electromagnetic fields, by the use of a combination of the Lorentz force law and the Liénard-Wiechert potential. In particular, one can substitute the explicit expressions for the $\mathbf{E}$ and $\mathbf{B}$ fields obtained from the Liénard-Wiechert formulas into the Lorentz force formula to derive the force between two charged particles moving on arbitrary paths in space. One can then derive a classical equation of motion for the particles using Newtonian mechanics, or its special relativistic correction.
Explicitly, we obtain this system of two ODEs, where $m_1, m_2, q_1, q_2$ are the masses and charges of the two particles, $\mathbf{r}_1(t), \mathbf{r}_2(t)$ are the paths, and other quantities are defined in the Wikipedia page for the Liénard-Wiechert potential:
$$m_1{\mathbf{\ddot{r}}_1}(t) = {\frac {\mu_0c^2q_1q_2}{4\pi}}\left(1 + \left[\mathbf{\dot{r}}_1(t)\times\left[\frac{\mathbf{n}_2(t_r)}{c}\times\right]\right]\right)$$
$${\left[\frac{\mathbf {n}_2(t_r) -{\boldsymbol {\beta}_2(t_r)}}{\gamma_2 ^{2}(t_r)(1-\mathbf {n}_2(t_r) \cdot {\boldsymbol {\beta }_2(t_r)})^{3}|\mathbf {r}_1(t) -\mathbf {r}_2(t_r)|^{2}}+\frac{\mathbf {n}_2(t_r) \times {\big (}(\mathbf {n}_2(t_r) -{\boldsymbol {\beta }_2(t_r)})\times {{\boldsymbol {\dot{\beta} }_2(t_r)}}{\big )}}{c(1-\mathbf {n}_2(t_r) \cdot {\boldsymbol {\beta }_2(t_r)})^{3}|\mathbf {r}_1(t) -\mathbf {r}_2(t_r)|}\right]}$$
$$m_2{\mathbf{\ddot{r}}_2}(t) = {\frac {\mu_0c^2q_1q_2}{4\pi}}\left(1 + \left[\mathbf{\dot{r}}_2(t)\times\left[\frac{\mathbf{n}_1(t_r)}{c}\times\right]\right]\right)$$
$${\left[\frac{\mathbf {n}_1(t_r) -{\boldsymbol {\beta}_1(t_r)}}{\gamma_1 ^{2}(t_r)(1-\mathbf {n}_1(t_r) \cdot {\boldsymbol {\beta }_1(t_r)})^{3}|\mathbf {r}_2(t) -\mathbf {r}_1(t_r)|^{2}}+\frac{\mathbf {n}_1(t_r) \times {\big (}(\mathbf {n}_1(t_r) -{\boldsymbol {\beta }_1(t_r)})\times {{\boldsymbol {\dot{\beta} }_1(t_r)}}{\big )}}{c(1-\mathbf {n}_1(t_r) \cdot {\boldsymbol {\beta }_1(t_r)})^{3}|\mathbf {r}_2(t) -\mathbf {r}_1(t_r)|}\right]}$$
where $\mathbf {n}_2(t_r) = \frac{\mathbf{r}_1(t)-\mathbf{r}_2(t_r)}{|\mathbf{r}_1(t)-\mathbf{r}_2(t_r)|}$, $\boldsymbol {\beta}_2(t_r) = \frac{\mathbf{\dot{r}}_2(t_r)}{c}$ and $\gamma_2(t_r) = \frac{1}{\sqrt{1-|\boldsymbol {\beta}_2(t_r)|^2}}$, and similarly for $\mathbf {n}_1(t_r)$, $\boldsymbol {\beta}_1(t_r)$ and $\gamma_1(t_r)$. The retarded time is defined implicitly by the equation $t_r = t - \frac{1}{c}|\mathbf{r}_1(t)-\mathbf{r}_2(t_r)|$. I slightly abused notation by "factoring out" the cross product terms to avoid duplicating things, hopefully this is clear enough.
This can be generalized to a system of $N$ charged particles in a similar way. I have not performed the calculations myself due to the apparent complexity of the resulting equation of motion, but in principle one could check whether the solutions correspond to what is predicted with the electron spiraling down into the nucleus, or whether it leads to something close to elliptical orbits.
My intuition tells me that we won't observe the kind of spiralling down predicted by assuming energy is stored in the electric fields throughout all of space in this setup. Instead we consider the $\mathbf{E}$ and $\mathbf{B}$ fields as useful mathematical abstractions to simplify the expression given above into more manageable components. The interpretation of the Poynting vector would be as energy flux density which would exist only if there exist other charges that would get accelerated by the fields. In particular, the atomic stability problem would require additional charges being present near the hydrogen atom, which would obviously affect the electron's orbit through additional force terms as a multi-body problem. In that scenario, there would be an energy transfer between the electron and other charges nearby. But even then it is not clear that the electron automatically loses energy, because the nearby particles will in turn radiate and couple with the electron's motion.
We can simplify the equations above by assuming the mass of particle 2 to be very large, so that it remains effectively stationary in some inertial frame of reference, and located at the origin. Then the equations above simplify greatly and we have the following equation of motion for particle 1 (with $\mathbf{r}_2(t) = \mathbf{0}$):
$$m_1{\mathbf{\ddot{r}}_1}(t) = {\frac {\mu_0c^2q_1q_2}{4\pi}}\frac{\mathbf {r}_1(t)}{|\mathbf {r}_1(t)|^{3}}$$
which is just the equation of motion of a charged particle moving in an electrostatic potential with the Coulomb force. This limiting model is formally identical to the model of a planet orbiting a massive object under Newtonian gravitation, and we clearly have elliptical orbits. Therefore it would be incorrect to claim that classical electromagnetism predicts instabilities of the atom (in this limiting case with very massive nucleus, at least) if we don't consider the energy supposedly stored in the EM fields. Also it deals away with the self-energy problem of a charged particle (integrating the "electrical field energy density" over all spaces gives an infinite result, which would be simply a meaningless calculation since there is no actual energy stored in such a field).
I hope what I said above was sufficiently clear, and I would be curious whether solutions to the equations of motion I have described have been calculated (or approximated somewhat) to predict classical orbits of an electron around the nucleus.
Note also that I am not questioning the validity of quantum mechanics and more detailed theories of matter. I am simply wondering whether the specific problem of atomic instability supposedly predicted by classical electromagnetics arises only due to the assumed existence of electromagnetic fields carrying energy, or whether a field-free formulation of electromagnetics using the equations of motion above is also subject to this problem. I am sure there are other problems that this model cannot resolve, such as the existence of discrete atomic emission and absorption spectra. But the important observation I wanted to make is that starting from the classical Maxwell equations and Lorentz force, we can derive the Liénard-Wiechert potential, and then derive the explicit equations of motion above, and finally forget about the existence of the $\mathbf{E}$ and $\mathbf{B}$ fields. This leads to a classical model of a two-body atom with stable orbits.
 A: The essential conceptual problem in your treatment is the fact that you assume that the radiation should somehow emerge solely from particle $A$ acting on particle $B$. This is not true, the radiation (and radiation-reaction) comes from the particle $B$ acting on itself! To understand this, you must consider $A$ and $B$ to be finite bodies first. Then the radiation schematically emerges as:


*

*all the charges in the body $B$ "sending out" their own electromagnetic potential on their null cones $|\vec{r}-\vec{r}_B|-ct = 0$, 

*body $A$ accelerating charges within body $B$,

*and finally, the accelerated charges within $B$ interacting with some of the electromagnetic potential from within $B$ (point 1.)!
(You can swap $B$ and $A$ to get the radiation from particle $A$ as well.) 
When the dust settles, you can take a limit of the sizes of the bodies going to zero and you get the famous Abraham-Lorentz-Dirac force acting on either $A$ or $B$:
$$F^{\mathrm{ALD}}_\mu = \frac{\mu_o q^2}{6 \pi m c}
\left[\frac{d^2 p_\mu}{d \tau^2}-\frac{p_\mu}{m^2 c^2}
\left(\frac{d p_\nu}{d \tau}\frac{d p^\nu}{d \tau}\right)
\right]$$
Nevertheless, the derivation of this result is notoriously challenging both conceptually and technically. The reason for that is that when you treat the body $B$ as an infinitely small particle, it should not be able to interact with data on its own light-cone because that would mean it is moving beyond the speed of light! On the other hand, the particle is on its own lightcone at $t=0$ and $\vec{r} = \vec{r}_B$, and the potential and the Lorentz force diverge at that exact position!
The only way to resolve this in a rigorous way is to assume, as already mentioned above, that the bodies in questions are of finite spatial extent and finite charge densities. Then, you take the limit of the size going to zero. This was carefully redone in 2009 by Gralla, Harte & Wald and I recommend that paper for further information. (The reason why this topic has received heightened interest recently is the fact that gravitational-wave inspirals of small stellar-mass astrophysical objects into supermassive black holes can be treated exactly in the approximation of a "self-forced particle", see Barack & Pound, 2018.)
You can derive the Larmor formula from a particular perturbative approximation of the ALD force called order reduction. First you take the $\mathrm{d}p^\mu/\mathrm{d}\tau$ for the particle without radiation reaction and insert it into $F^{\rm ALD}_\mu$. Then the Larmor formula is just the rate with which this force is taking energy from the particle. 

EDIT: A broader historical discussion
Jan Lálinský reminded me that there exist formulations of classical electrodynamics that 1) agree with most of the predictions of Maxwell equations in the continuum limit (given a certain "absorbing universe" postulate), and 2) where the "point particle" is not the limit of a finite body and does not feel any self-force. A brief review of these "Schwarzschild-Tetrode-Fokker(-Frenkel)" (STF) electrodynamics was given by Wheeler & Feynman in 1949.
Depending on how exactly do you implement the "absorber universe", the quasi-neutral ensemble of particles far away from your system, the planetary atom is also usually unstable in STF electrodynamics. This is because the energy tends to be stolen from the atom by the ensemble of far-away particles and dissipated (even though possibly at a slower rate). On one hand, this is a nice "Machian" twist on electrodynamics, since the notion of the field is emergent from the physical particles, and the particles would not radiate energy if there were no other particles to pass the energy to. On the other hand, the STF electrodynamics tend to have curious non-local properties such as the absorber universe "knowing" about an action on the particle an infinite time before the action itself occurs! This makes the theory physically unsatisfactory to me.
Consider the following example of a pulsar, whose pulse we detect and whose rotation rate slows down as a consequence of the radiative energy loss. In mainstream electrodynamics, we talk about electromagnetic waves traveling for eons through space from the pulsar as independent energy-carrying entities, while the STF theory defies this picture. While in the mainstream electrodynamics the wave took the energy away from the pulsar and made it slow down its rotation rate, in the STF theory the pulsar slows down (or not) thanks to the fact that it "knows" that energy-receiving objects such as your antenna will be there in a thousand years!!!
Ultimately, both the usual particle+field and STF theories are wrong, and the correct theory of electrodynamics of fundamental point particles is quantum electrodynamics (and even more ultimately the Standard model), so this is more of an academic discussion. However, I find the STF picture grossly undidactic as compared to the understanding of classical electrodynamics as the theory of the electromagnetic field sourced by finite continua that we sometimes limit towards approximate point particles.
A: John Lighton Synge had similar idea and he analyzed numerically the equations of motions for two oppositely charged particles of arbitrary masses where only retarded EM forces are present.
J. L. Synge, On the electromagnetic two–body problem., Proc. Roy. Soc. A 177 118–39 (1940)
https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1940.0114
He found that the system still collapses, but the greater the difference in masses, the more slowly it does so. For hydrogen atom, the time of collapse turned out to be hundreds of times longer than the time naively obtained from the Larmor formula.
The collapse is due to simple but perhaps too simple assumption: that the force acting on any particle is just retarded EM force due to the other particle. Then, because the motion of the particles is anticorrelated (when one moves left, the other moves right), the system radiates EM energy away.
If we introduce into the model background EM radiation that acts on both particles (additional forces), collapse ceases to be inevitable, because the radiated energy can be supplied by the background radiation forces. There are some papers on that - for more on this, see also my answer here:
The classical electrodynamic atom
A: This subject has long been an interest of mine and from what I can find, it seems to me that the theory of electromagnetic direct particle interaction is not well-developed enough to answer your question. The reason is that the equations you're working with are not just any kind of ODEs, they are delay differential equations. As far as I know, there is no general solution of 2 particles directly interacting in this way for 2 dimensions or more. I think in 1 dimension the problem only has a global solution if the charges have the same sign (though I could be mistaken and the attractive case might be solved too).
There is one group I've found tackling this question in an interesting way. They work with a formalism where they assume direct particle interaction along the light cone AND Maxwell's equations for the field (in which the fields are uniquely prescribed by the particles' trajectories and not independent, dynamic degrees of freedom). They have some results showing that some solutions to this system are also solutions to the direct particle interaction system. Their approach is very mathematical but I will provide links to their work here for your interest, though I am not qualified to vet their approach:
https://iopscience.iop.org/article/10.1088/1751-8113/49/44/445202/pdf
https://www.tandfonline.com/doi/abs/10.1080/03605302.2013.814142
https://www.sciencedirect.com/science/article/pii/S0022039616000243
https://arxiv.org/abs/1603.05115
From my read of things: If you consider a universe of only two particles, their orbits would be completely stable. Any instability comes from the interaction of the two particle system with a larger bath of many, many particles. This phenomenon can be captured by the Dirac-Abraham-Lorentz force, which arises from these complicated many-particle interactions. In the standard theory, this can be added on as a sort of fudge factor to Maxwell's equations, but in doing so a ton of complications are introduced and mathematically the resulting ODEs may not be well-posed. Nonetheless, the standard Maxwell-Lorentz theory forms the basis of quantization that results in QED, but one can quantize direct particle theories instead. The result is a different quantum field theory of the electromagnetic field with its own pros and cons (the cons outweigh the pros I believe or else it would be the standard approach).
A: I agree with Void's answer, but I'll offer another side: By taking the mass of one particle to be much larger than the other, you've ended up with one charge doing all the radiating under the influence of the static Coulomb field of the other: your model of an electron orbiting a nucleus has the same equation of motion as a small body orbiting a much large body in a Newtonian gravitational system. Your model doesn't predict spiralling into the nucleus because you've used the standard Lorentz force without the radiation damping term. It's mathematically correct, but violates conservation of energy and momentum of the whole system when radiation is significant. 
Dirac$^1$ addressed this problem by deriving an equation of motion for an arbitrarily moving charge using the local conservation of energy and momentum for a tube surrounding the charge:
$$1/2q^2\epsilon^{-1}\dot{v}_{\mu} - qv_{\nu}f^{\nu}_{\mu} = \dot{B}_{\mu}$$
Where $q$ is the charge, $\epsilon$ the radius of the tube, $v$ the four-velocity, $f$ the field bounded to the charge, $B$ an undetermined four-vector.
To get further, he had to make further assumptions on how simple the equation was likely to be, and add a negative mass to compensate the Coulomb electromagnetic mass contribution $\rightarrow \infty$ as $\epsilon \rightarrow 0$, getting:
$$m\dot{v}_{\mu} - 2/3q^2\ddot{v}_{\mu} - 2/3q^2\dot{v}^2 {v}_{\mu} = ev_{\nu}F^{\nu}_{\mu\;\text{in}}$$
Dirac's derivation has the advantage of ignoring the structure of the charge and the Poincare stresses holding it together; unlike earlier models created by Lorentz, Abraham and Schott. However, it has problems with pre-acceleration and causality leading to it being modified by the Landau-Lifshitz equation.

[1] Dirac, P.A.M. Proc. R. Soc. London A 167, 148 (1938).
A: 
The full equations are complicated so I take a limit as the mass of the nucleus goes to infinity, and (assuming my calculations are correct), I obtain a simpler equation of motion formally identical to that of planetary motion around the sun. This predicts elliptical trajectories for most initial conditions.

From a comment of the OP:

The question is whether this classical model predicts the instability of the atom and whether it is truly a consequence of classical electromagnetism

The following discusses the stability part:
The basic question is whether the solutions you find are stable, or metastable, i.e. a small perturbation, like, radiation in the field of the other, or  atomic vibrations, will send the electron down to the nucleus. 
(I do remember that metastable states can exist in the classical solutions, but cannot find the reference)
From glancing through your derivation , I cannot understand what you do with radiation, i.e. how you could perturb your solution to see if it is stable or metastable. Radiation is an experimental fact. A charge radiates energy away in the electromagnetic spectrum when accelerating in a field . This is an experimental fact. Where is radiation in your formulas?
Bohr did obtain planetary solutions, but needed to impose quantization of angular momentum in order to have stability. ( radiation would carry of a unit of angular momentum)
I suspect that this is the case with your solutions, they are metastable not taking into account radiation.
