Because this derivation is often given as homework in classes that teach relativity, I will not show step-by-step derivations. Instead, I will only show the steps that books tend to give. If you are confident in your use of tensor notation, you shouldn't have a problem filling in the steps. However, it can get extremely messy when it calls for a change in indices. I will again reference "A General Relativity Workbook" by Thomas A. Moore.
We use tensors in this derivation so that it is valid in arbitrary coordinates. Let's go back to our particles and define their positions by $ x^{\alpha}(\tau) $ and $ \bar{x}^{\alpha}(\tau) \equiv x^{\alpha}(\tau) + n^{\alpha}(\tau) $, where $ \tau $ is our proper time and $ \textbf{n} $ is our infinitesimal separation four-vector. Keep in mind that although this derivation is similar to that in Newtonian mechanics, this derivation uses relativistic concepts. Just as with the Newtonian derivation, we must write expressions that govern the motion of both particles along their geodesics.
$$ 0 = \frac{d^2 x^{\alpha}}{d \tau^2} + \Gamma^{\alpha}_{\mu \nu} \frac{dx^{\mu}}{d \tau} \frac{dx^{\nu}}{d \tau} , \qquad 0 = \frac{d^2 \bar{x}^{\alpha}}{d \tau^2} + \bar{\Gamma}^{\alpha}_{\mu \nu} \frac{d \bar{x}^{\mu}}{d \tau} \frac{d \bar{x}^{\nu}}{d \tau} $$
Here, we express the particles' coordinate accelerations in terms of the Christoffel symbols instead of having them in terms of the gravitational potential. Christoffel symbols are combinations of first derivatives of the metric that describe effects of parallel transport in manifolds. Just like we did in the Newtonian derivation, we can use Taylor series to expand Christoffel symbol of the second particle around the value of the first particle since the separations are infinitesimal (shown below).
$$ \bar{\Gamma}^{\alpha}_{\mu \nu} (at \: \bar{x}^{\alpha} (\tau)) \approx \Gamma^{\alpha}_{\mu \nu} (at \: x^{\alpha} (\tau)) + n^{\sigma} [ \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu}] (at \: x^{\alpha} (\tau)) $$
Now, we know $ \bar{x}^{\alpha}(\tau) \equiv x^{\alpha}(\tau) + n^{\alpha}(\tau) $ so we will plug this and the expansion above into the equation of motion of the second particle which gives us
$$ 0 = \frac{d^2 (x^{\alpha} + n^{\alpha})}{d \tau^2} + [\Gamma^{\alpha}_{\mu \nu} + n^{\sigma} ( \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu})] \frac{d (x^{\mu} + n^{\mu})}{d \tau} \frac{d (x^{\nu} + n^{\nu})}{d \tau} $$
$$ \qquad \qquad \; \; = \frac{d^2 x^{\alpha}}{d \tau^2} + \frac{d^2 n^{\alpha}}{d \tau^2} + [\Gamma^{\alpha}_{\mu \nu} + n^{\sigma} ( \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu})] \left( \frac{d x^{\mu}}{d \tau} + \frac{d n^{\mu}}{d \tau} \right) \left( \frac{d x^{\nu}}{d \tau} + \frac{d n^{\nu}}{d \tau} \right) $$
We want to simplify this more by multiplying out everything and using our geodesic equations. Doing this will give us
$$ ** \quad 0 = \frac{d^2 n^{\alpha}}{d \tau^2} + 2 \Gamma^{\alpha}_{\mu \nu} u^{\mu} \frac{d n^{\nu}}{d \tau} + n^{\sigma} ( \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu} u^{\mu}) u^{\mu} u^{\nu} $$
where $ u^{\mu} \equiv \frac{d x^{\mu}}{d \tau} $ is the first particle's four-velocity. Because this is a four-vector we know that its derivative with respect to $ \tau $ is
$$ \left( \frac{d \textbf{n}}{d \tau} \right)^{\alpha} = \frac{d n^{\alpha}}{d \tau} + \Gamma^{\alpha}_{\mu \nu} u^{\mu} n^{\nu} $$
and so we take the second derivative and find it to be
$$ \left( \frac{d^2 \textbf{n}}{d \tau^2} \right)^{\alpha} = \left( \frac{d}{d \tau} \left[ \frac{d \textbf{n}}{d \tau} \right] \right)^{\alpha} = \frac{d}{d \tau} \left( \frac{d n^{\alpha}}{d \tau} + \Gamma^{\alpha}_{\mu \nu} u^{\mu} n^{\nu} \right) + \Gamma^{\alpha}_{\sigma \nu} u^{\sigma} \left( \frac{d n^{\nu}}{d \tau} + \Gamma^{\nu}_{\beta \gamma} u^{\beta} n^{\gamma} \right) $$
$$ = \frac{d^2 n^{\alpha}}{d \tau^2} + \frac{d \Gamma^{\alpha}_{\mu \nu}}{d \tau} u^{\mu} n^{\nu} + \Gamma^{\alpha}_{\mu \nu} \frac{d u^{\mu}}{d \tau} n^{\nu} + \Gamma^{\alpha}_{\mu \nu} u^{\mu} \frac{d n^{\nu}}{d \tau} + \Gamma^{\alpha}_{\sigma \nu} u^{\sigma} \frac{d n^{\nu}}{d \tau} + \Gamma^{\alpha}_{\sigma \nu} \Gamma^{\nu}_{\beta \gamma} u^{\sigma} u^{\beta} n^{\gamma} $$
Now, this is getting really messy so let us simplify by noticing two key things:
$$ 1. \frac{d \Gamma^{\alpha}_{\mu \nu}}{d \tau} = \frac{d x^{\sigma}}{d \tau} \frac{\partial \Gamma^{\alpha}_{\mu \nu}}{\partial x^{\sigma}} = u^{\sigma} \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu} $$
$$ 2. \; Changing \; \sigma \; to \; \mu \; will \; make \; two \; of \; the \; terms \; the \; same. $$
Using both of these in the messy equation we have so far simplifies it to
$$ \left( \frac{d^2 \textbf{n}}{d \tau^2} \right)^{\alpha} = \frac{d^2 n^{\alpha}}{d \tau^2} + \left( \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu} \right) u^{\sigma} u^{\mu} n^{\nu} + \Gamma^{\alpha}_{\mu \nu} \frac{d u^{\mu}}{d \tau} n^{\nu} + 2 \Gamma^{\alpha}_{\mu \nu} u^{\mu} \frac{d n^{\nu}}{d \tau} + \Gamma^{\alpha}_{\sigma \nu} \Gamma^{\nu}_{\beta \gamma} u^{\sigma} u^{\beta} n^{\gamma} $$
Ok, so now the tricky part comes. We need to do three things: use the geodesic equation of the first particle to eliminate $ \frac{d u^{\mu}}{d \tau} $, use the equation labeled with ** in order to eliminate $ \frac{d^2 n^{\alpha}}{d \tau^2} $, and rename indices to pull out a common factor of $ u^{\sigma} u^{\mu} n^{\nu} $. This will take a while. The trickiest part is definitely changing the indices because there are so many of them. Try to see what you can cancel out and remember the symmetries of the Christoffel symbols. After this, we will get
$$ \left( \frac{d^2 \textbf{n}}{d \tau^2} \right)^{\alpha} = \left( \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu} - \partial_{\nu} \Gamma^{\alpha}_{\mu \sigma} + \Gamma^{\alpha}_{\sigma \gamma} \Gamma^{\gamma}_{\mu \nu} - \Gamma^{\alpha}_{\nu \gamma} \Gamma^{\gamma}_{\mu \sigma} \right) u^{\sigma} u^{\mu} n^{\nu} $$
Since the left side is a tensor and everything outside the parenthesis on the right side is a tensor, then the stuff inside the parenthesis must be a tensor. Not surprisingly, the stuff inside is actually the Riemann tensor, but with the signs flipped. So the Riemann tensor, which is the measure of spacetime curvature, is defined as
$$ R^{\alpha}_{\mu \nu \sigma} \equiv \partial_{\nu} \Gamma^{\alpha}_{\mu \sigma} - \partial_{\sigma} \Gamma^{\alpha}_{\mu \nu} + \Gamma^{\alpha}_{\nu \gamma} \Gamma^{\gamma}_{\mu \sigma} - \Gamma^{\alpha}_{\sigma \gamma} \Gamma^{\gamma}_{\mu \nu} $$
Therefore, we can rewrite the last equation using the Riemann definition and we get
$$ \left( \frac{d^2 \textbf{n}}{d \tau^2} \right)^{\alpha} = - R^{\alpha}_{\mu \nu \sigma} u^{\sigma} u^{\mu} n^{\nu} $$
This is the Geodesic Deviation Equation. This is extremely useful because it allows us to see if spacetime is flat or curved, which we can't do by looking at the metric. So if the any of the components of the Riemann tensor are non-zero, then spacetime is curved.
