In addition to dmckee's answer this link summarizes the experiments of antiprotons catching protons and neutrons, creating temporary nuclei.
It is the symmetric state to the one in the question : an anti-proton-neutron nucleus that lasts for a bit (fig 5.7). Anti-protons can be made in abundance and controlled experimentally because they are charged, anti-neutrons do not give easy handles.
The relativistic generalization of the formula; and the replacement of the electron mass by the reduced mass are clearly two basically independent steps (at least in the leading approximation), and both of them have to be applied to agree with the experiment.
The most accurate experimental value of the energy is $-13.59844\,{\rm eV}$, see NIST:
http://webbook.nist.gov/cgi/cbook.cgi?ID=C12385136&Units=CAL&Mask=20
This takes the relativistic corrections into account, as well as the motion of the nucleus (i.e. the replacement of the electron mass by the reduced mass). Your value $-13.60587\,{\rm eV}$ was a theoretical value calculated with the electron mass, not reduced mass (although someone may have used slightly different values of the fine-structure constant etc. than you have, therefore the tiny discrepancy), so this value can't be expected to match the experiment too accurately.
If one goes to higher orders, these two effects (relativistic; proton's motion) cease to be truly independent of each other, and higher-order terms generally exist that can't be incorporated by a simple replacement of $m_e$ by $\mu$ or by the simple relativistic square roots.
The exact calculations at that precision actually require the full quantum electrodynamics, see e.g.
https://arxiv.org/abs/quant-ph/0511197
but only the most precise experiments may be sensitive to these corrections. The binding energy isn't measured too accurately. On the other hand, the transition energies are measured much more accurately due to spectroscopy and they're sufficiently good checks to verify whether one has incorporated all the corrections – and their mutual interactions – correctly.
Note that in the exact treatment, there are additional small terms that have to be taken into account to get "truly precision values" to match the experiment, such as the interaction of the electron's and nucleus' magnetic moment; and the nonzero size of the nucleus that truncates the Coulomb potential for small enough $r$ (distance of the electron and the nucleus).
Best Answer
Consider the ratio of their shifts in energy. Since all the numbers are the same except for the reduced mass you are looking at something like: $$ \frac{(E_i -E_f)_p}{(E_i -E_f)_d}=\frac{\mu_p}{\mu_d}=\frac{m_p m_e(m_d +m_e)}{(m_p +m_e)m_d m_e} $$ If we cancel the $m_e$ and then pull a $m_d$ out of the top and a $m_p$ out of the bottom we get $$ \frac{(E_i -E_f)_p}{(E_i -E_f)_d}=\frac{(1+\frac{m_e}{m_d})}{(1+\frac{m_e}{m_p})} $$ I hope this helps.