The aim of the isotropic coordinates is to write the metric in the form where the spacelike slices are as close as possible to Euclidean. That is, we try to write the metric in the form:
$$ ds^2 = -A^2(r)dt^2 + B^2(r)d\Sigma^2 $$
where $d\Sigma^2$ is the Euclidean metric:
$$ d\Sigma^2 = dr^2 + r^2(d\theta^2 + \sin^2\theta d\phi^2) $$
So let's use the substitution $r\rightarrow r'$ and write down our metric:
$$ ds^2 = -\left(1-\frac{2M}{r'}\right)dt^2 + B^2(r')\left(dr'^2 + r'^2(d\theta^2 + \sin^2\theta d\phi^2)\right) $$
If we compare this with the Schwarzschild metric:
$$ ds^2 = -\left(1-\frac{2M}{r}\right)dt^2 + \frac{dr^2}{1 - 2M/r} + r^2(d\theta^2 + \sin^2\theta d\phi^2) $$
Then for the angular parts to be equal we must have:
$$ B^2(r')r'^2 = r^2 $$
And for the radial parts to be equal we must have:
$$ B^2(r')dr'^2 = \frac{dr^2}{1 - 2M/r} $$
Divide the second equation by the first to eliminate $B$ and we end up with:
$$ \frac{dr'^2}{r'^2} = \frac{dr^2}{r^2 - 2Mr} $$
And then just take the square root and integrate and we get the substitution you describe:
$$ r = r'\left(1 + \frac{M}{2r'}\right)^2 $$
There is a relatively fast approach to computing the Riemann tensor, Ricci tensor and Ricci scalar given a metric tensor known as the Cartan method or method of moving frames. Given a line element,
$$ds^2 = g_{\mu\nu}dx^\mu dx^\nu$$
you pick an orthonormal basis $e^a = e^a_\mu dx^\mu$ such that $ds^2 = \eta_{ab}e^a e^b$. The first Cartan structure equation,
$$de^a + \omega^a_b \wedge e^b = 0$$
allows one to solve for the spin connection components $\omega^a_b$ from which one can compute the Ricci tensor in the orthonormal basis:
$$R^a_b = d\omega^a_b + \omega^a_c \wedge \omega^c_b.$$
The entire process simply requires exterior differentiation of the basis and spin connection. The Riemann components may be deduced from the relation,
$$R^a_b = R^a_{bcd} \, e^c \wedge e^d$$
possibly with a factor of $\frac12$ depending on your conventions. To convert back to the coordinate basis, one must simply contract with the basis back:
$$R^\mu_{\nu \lambda \kappa} = (e^{-1})^\mu_a \, R^a_{bcd}\, e^b_\nu \, e^c_\lambda \, e^d_\kappa.$$
For an explicit calculation see my previous answers here, here and here. The gravitational physics lectures at pirsa.org also provide explicit examples. As for using computer algebra systems, if all you're looking to do is compute curvature tensors, Hartle's textbook for Mathematica is your best option or the GREAT package. If you'd like to do more advanced stuff like perturbation theory, then xAct is required.
Best Answer
According to Mathematica, and assuming I haven't made any silly errors typing in the metric, I get the non-zero components of $R^\mu{}_{\nu\alpha\beta}$ to be: