The simplest way to explain the Christoffel symbol is to look at them in flat space. Normally, the laplacian of a scalar in three flat dimensions is:
$$\nabla^{a}\nabla_{a}\phi = \frac{\partial^{2}\phi}{\partial x^{2}}+\frac{\partial^{2}\phi}{\partial y^{2}}+\frac{\partial^{2}\phi}{\partial z^{2}}$$
But, that isn't the case if I switch from the $(x,y,z)$ coordinate system to cylindrical coordinates $(r,\theta,z)$. Now, the laplacian becomes:
$$\nabla^{a}\nabla_{a}\phi=\frac{\partial^{2}\phi}{\partial r^{2}}+\frac{1}{r^{2}}\left(\frac{\partial^{2}\phi}{\partial \theta^{2}}\right)+\frac{\partial^{2}\phi}{\partial z^{2}}-\frac{1}{r}\left(\frac{\partial\phi}{\partial r}\right)$$
The most important thing to note is the last term above--you now have not only second derivatives of $\phi$, but you also now have a term involving a first derivative of $\phi$. This is precisely what a Christoffel symbol does. In general, the Laplacian operator is:
$$\nabla_{a}\nabla^{a}\phi = g^{ab}\partial_{a}\partial_{b}\phi - g^{ab}\Gamma_{ab}{}^{c}\partial_{c}\phi$$
In the case of cylindrical coordinates, what the extra term does is encode the fact that the coordinate system isn't homogenous into the derivative operator--surfaces at constant $r$ are much larger far from the origin than they are close to the origin. In the case of a curved space(time), what the Christoffel symbols do is explain the inhomogenities/curvature/whatever of the space(time) itself.
As far as the curvature tensors--they are contractions of each other. The Riemann tensor is simply an anticommutator of derivative operators--$R_{abc}{}^{d}\omega_{d} \equiv \nabla_{a}\nabla_{b}\omega_{c} - \nabla_{b}\nabla_{a} \omega_{c}$. It measures how parallel translation of a vector/one-form differs if you go in direction 1 and then direction 2 or in the opposite order. The Riemann tensor is an unwieldy thing to work with, however, having four indices. It turns out that it is antisymmetric on the first two and last two indices, however, so there is in fact only a single contraction (contraction=multiply by the metric tensor and sum over all indices) one can make on it, $g^{ab}R_{acbd}=R_{cd}$, and this defines the Ricci tensor. The Ricci scalar is just a further contraction of this, $R=g^{ab}R_{ab}$.
Now, due to Special Relativity, Einstein already knew that matter had to be represented by a two-index tensor that combined the pressures, currents, and densities of the matter distribution. This matter distribution, if physically meaningful, should also satisfy a continuity equation: $\nabla_{a}T^{ab}=0$, which basically says that matter is neither created nor destroyed in the distribution, and that the time rate of change in a current is the gradient of pressure. When Einstein was writing his field equations down, he wanted some quantity created from the metric tensor that also satisfied this (call it $G^{ab}$) to set equal to $T^{ab}$. But this means that $\nabla_{a}G^{ab} =0$. It turns out that there is only one such combination of terms involving first and second derivatives of the metric tensor: $R_{ab} - \frac{1}{2}Rg_{ab} + \Lambda g_{ab}$, where $\Lambda$ is an arbitrary constant. So, this is what Einstein picked for his field equation.
Now, $R_{ab}$ has the same number of indicies as the stress-energy tensor. So, a hand-wavey way of looking at what $R_{ab}$ means is to say that it tells you the "part of the curvature" that derives from the presence of matter. Where does this leave the remaining components of $R_{abc}{}^{d}$ on which $R_{ab}$ does not depend? Well, the simplest way (not COMPLETELY correct, but simplest) is to call these the parts of the curvature derived from the dynamics of the gravitational field itself--an empty spacetime containing only gravitational radiation, for example, will satisfy $R_{ab}=0$ but will also have $R_{abc}{}^{d}\neq 0$. Same for a spacetime containing only a black hole. These extra components of $R_{abc}{}^{d}$ give you the information about the gravitational dynamics of the spacetime, independent of what matter the spacetime contains.
This is getting long, so I'll leave this at that.
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
The idea is that we want to define some notion of curvature for a manifold that intuitively agrees with the intuition we have about curvature.
The genius insight that leads to the desired definition is the notion of parallel transport. Speaking non-rigorously here, the basic idea is that if you transport a tangent vector on a manifold parallel to itself all the way around a closed curve, then the vector will come back to itself in flat spaces, but it will become a different vector in a curved space.
To see why the notion of parallel transport has anything to do with curvature, think for example, of the Euclidean plane $\mathbb R^2$ versus the two-dimensional sphere $S^2$.
Consider the curve consisting of an equilateral triangle with one vertex at the origin. Now imagine placing a vector emanating from the origin, and imagine moving that vector along the triangle, keeping its "tail" on the triangle, and making sure to keep the vector parallel to itself the whole time. If you transport the vector once around the triangle back to the origin in this way, then you get the same vector back.
Something drastically different happens if you do the same thing on the sphere as the following diagram from the wiki page on parallel transport indicates
If you move a vector from point A back to itself along the curve indicated in the diagram, the vector does not return to itself. This happens because the sphere is curved.
In fact, the notion of parallel transport can be used to completely characterize what we mean by curvature. The logic you'll find in many books on GR and differential geometry is roughly as follows:
Define the notion of a connection (basically this defines what you mean by taking derivatives on the manifold).
Use the connection to define the notion of parallel transport which agrees with our intuition of parallel transport in, for example, the sphere example above.
Show that there is a tensor that measures precisely how much the components of a vector change when it is parallel transported along a small closed curve on the manifold.
Call this tensor the Riemann tensor, and use it as the object that captures the notion of curvature.
There is a great discussion of this in a lot of books. I personally like the discussion on pages 36-38 of Wald's General Relativity.
Addendum. Wald actually shows that if you consider a curve bounding a small two-dimensional patch parameterized by coordinates $s$ and $t$ on the given manifold, then the change $\delta v^a$ in the components of a vector transported along the boundary of this patch satisfies \begin{align} \delta v^a = \delta t\,\delta s\, v^dT^cS^bR_{cbd}^{\phantom{cbd}a} \end{align} where $\delta t\,\delta s$ is the area of the patch, and $T^c$ and $S^b$ are the tangents to the curves of constant $s$ and $t$ respectively.