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).
While you should learn to compute the integrals, let me point out in this additional answer that one can compute this without evaluating any integrals.
The Hamiltonian is:
$$H = \frac{p^2}{2m} + \frac{e^2}{r}$$
the energy eigenvalue of the $\left|n,l,m\right>$ eigenstate is
$$E_{nlm} = \frac{m e^4}{2\hbar^2 n^2}\tag{1}$$
here I've written this explicitly in terms of the constants that appear in the Hamiltonian. If we add the term $$V = \frac{g}{r}$$ as a perturbation to the Hamiltonian, then we know from first order perturbation theory that the perturbed energy eigenvalues are given to first order in $g$ by:
$$E(g) = E_{nlm} + g \left<nlm\right|\frac{1}{r}\left|nlm\right>+\mathcal{O}(g^2)\tag{2}$$
We can also compute the perturbed energy directly by modifying the charge $e$. We can write:
$$H + V = \frac{p^2}{2m} + \frac{e'^2}{r}$$
with:
$$e'^2 = e^2 + g$$
From Eq. (1) it follows that:
$$E(g) = \frac{me'^4}{2\hbar^2} = E_{nlm} + g\frac{me^2}{\hbar^2 n^2} + \mathcal{O}(g^2)$$
Comparing this to Eq. (2) yields:
$$\left<nlm\right|\frac{1}{r}\left|nlm\right> = \frac{me^2}{\hbar^2 n^2}=\frac{1}{a_0 n^2}$$
Best Answer
Your problem lies in assuming that $$ \nabla^2 = \frac{\partial^2}{\partial r^2} + \cdots $$
This is not the case, you need to use $$ \nabla^2 = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\cdots $$ Then will you obtain the correct answer of $-1/2$.