I think there is a method that I believe is rather simple. Take a look:
There is a thing called 'normal Riemann coordinates'. In this coordinates the metric is expanded around the origin, and the coefficients of expansion are expressed in terms of the Riemann tensor. I suggest that you read about them and check whether the coordinates described below are normal. All needed information can be found here.
Take the coordinates $x_i,\,x=\sqrt{x_i x_i}$ to be (I write $r$ so it cant be confused with $R$ the contraction of Ricci tensor):
$$
\psi=x/r\\
x_1=r\psi\sin\theta\cos\phi\\
x_2=r\psi\sin\theta\sin\phi\\
x_3=r\psi\cos\theta
$$
Compare it with usual spherical coordinates in $\mathbb{R}^3$, we know:
$$
dx_idx_i=r^2d\psi^2+r^2\psi^2(d\theta^2+\sin^2\theta d\phi^2)
$$
So we can rewrite your metric as
$$
ds^2=\frac{\sin^2\psi}{\psi^2}dx_i dx_i+r^2d\psi^2\left(1-\frac{\sin^2\psi}{\psi^2}\right)=\\
=r^2\frac{\sin^2(x/r)}{x^2}dx_i dx_i+\frac{(x_idx_i)^2}{x^2}\left(1-r^2\frac{\sin^2(x/r)}{x^2}\right)
$$
Now, if we believe that $x_i$ are normal, then we have:
$$
g_{ik}(x)=\delta_{ik}-\frac{1}{3}R_{iakb}x^ax^b+...
$$
So you can just expand the metric and find the symmetrized form $\frac{1}{2}\left(R_{iakb}+R_{ibka}\right)$ of the Riemann tensor at the origin that is sufficient to contract it to the Ricci tensor. The good point is that it is the same level of complexity (in fact, the 'normal' metric of the same form) in any dimension. As soon as you will find $R_{ik}\sim g_{ik}$ in this coordinate system, it holds in any coordinate system. And, due to the symmetry of the sphere, it also holds at every point.
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
Here's a way to find the Riemann tensor of the 3-sphere with a lot of intelligence but no calculations.
At any point $p$ on a sphere, all directions look the same. Therefore there can be no privileged vector at a point $p$. Now consider the eigenvalue problem for the Ricci tensor, $$R^\alpha{}_\beta x^\beta = \lambda x^\alpha.$$ Since no vector is better than any other vector, either no vector is an eigenvector, or every vector is. Since the Ricci tensor is symmetric, it must be the latter. But the only way for every vector to be an eigenvector is if the tensor is proportional to the identity tensor. Hence, $$R^\alpha{}_\beta = \lambda \delta^\alpha{}_\beta \Rightarrow R_{\alpha\beta} = \lambda g_{\alpha\beta}.$$ Now take the trace of both sides to conclude that $\lambda = R/3$ where $R$ is the curvature scalar.
Now we need to be even more clever. A diagonal element $R_{\alpha\beta}x^\alpha x^\beta$ of the Ricci tensor in $n$ dimensions is $(n-1)$ times the average of the sectional curvature over planes containing $x^\alpha$, for $x_\alpha x^\alpha = 1$. The important thing is that the sectional curvature is the Gaussian curvature of the 2-surface generated by geodesics starting at $p$. What can this 2-surface be? It's a 2-surface characterized by constant curvature (because of all the symmetry we have) -- that's (some part of) a 2-sphere! The Gaussian curvature of a 2-sphere is $1/r^2$. Thus $2/r^2 = \lambda g_{\alpha\beta}x^\alpha x^\beta = \lambda$ and we can conclude $$R_{\alpha\beta} = \frac{2}{r^2} g_{\alpha\beta}.$$