As you have it set up right now, you should be able to take as many derivatives as you want and evaluate it at $t=0$ just by using your initial conditions which is quite brilliant. The only downside to Taylor expansions is that the approximation will quickly go bad as $t$ increases. You also typically get astonishingly diminishing returns as you add more terms to the series. Disclaimer is now over, so it's time to get on to your question.
Let's start by rewriting your ODE a little bit. Let's define a vector $\Delta \vec{r}_k$ so we can write
$$\begin{align} &\vec{r}_k(t)-\vec{r}_i(t) = \Delta \vec{r}_k \\\\\implies
& \|\vec{r}_k(t)-\vec{r}_i(t) \|^2 = \Delta \vec{r}_k \cdot \Delta \vec{r}_k\\\\\implies &
\frac{d^2}{dt^2}\vec{r_i}(t)=G
\sum_{k=1}^{n}
m_k(\vec{r}_k(t)-\vec{r}_i(t)) \:
\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2}
\end{align}$$
where $ \: \cdot \: $ is of course the dot product between the vectors. Now we can take the product and chain rules to the statement quite a bit more easily. Observe
$$ \begin{align} \frac{d^3}{dt^3}\vec{r_i}(t)&=\frac{d}{dt}\bigg[ G
\sum_{k=1}^{n}
m_k(\vec{r}_k(t)-\vec{r}_i(t)) \:
\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2}\bigg]\\\\
&= G
\sum_{k=1}^{n}
m_k\frac{d}{dt}\bigg[\vec{r}_k(t)-\vec{r}_i(t)\bigg] \:
\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2}\\
\end{align} $$
$$ + m_k(\vec{r}_k(t)-\vec{r}_i(t)) \:
\frac{d}{dt}\bigg[\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2} \bigg] \text{ .} $$
Here it is quite clear that
$$ \frac{d}{dt}\bigg[\vec{r}_k(t)-\vec{r}_i(t)\bigg] = \frac{d}{dt}\vec{r}_k(t)-\frac{d}{dt}\vec{r}_i(t) = \frac{d}{dt} \Delta \vec{r}_k$$
but it may not be clear given your background how to compute
$$ \frac{d}{dt}\bigg[\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2} \bigg] \text{ .} $$
Remember that the dot product between vectors produces a scalar so we are really just applying the chain rule to a scalar function. Observe
$$ \frac{d}{dt}\bigg[\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2} \bigg] = -\frac{3}{2}\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-5/2} \: \frac{d}{dt}\bigg[\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \bigg] \text{ .}$$
Using upper indices to denote the component of the vector, the dot product will just look something like
$$ \begin{align} & \Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) = (\Delta \vec{r}_k^{\:1})^2 + (\Delta \vec{r}_k^{\:2})^2 + \cdots \\\\\implies &
\frac{d}{dt}\bigg[\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \bigg] = 2\Delta \vec{r}_k^{\: 1}\frac{d}{dt}\bigg[\Delta \vec{r}_k^{\: 1}\bigg] + 2\Delta \vec{r}_k^{\: 2}\frac{d}{dt}\bigg[\Delta \vec{r}_k^{\: 2}\bigg] + \cdots \\\\\implies &
\frac{d}{dt}\bigg[\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \bigg] = 2\Delta \vec{r}_k \cdot \frac{d}{dt}\bigg[\Delta \vec{r}_k\bigg]\end{align} $$
so now we can answer the original question. The next derivative of your ODE is
$$ \begin{align} \frac{d^3}{dt^3}\vec{r_i}(t)&= G
\sum_{k=1}^{n}
m_k\frac{d}{dt}\bigg[\Delta\vec{r}_k(t)\bigg] \:
\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-3/2}\\
\end{align} $$
$$ + m_k\Delta\vec{r}_k(t) \:
\bigg( -\frac{3}{2}\left(\Delta\vec{r}_k(t) \cdot \Delta\vec{r}_k(t) \right)^{-5/2} \: \bigg)\bigg(2\Delta \vec{r}_k \cdot \frac{d}{dt}\bigg[\Delta \vec{r}_k\bigg]\bigg) $$
which honestly is not as bad as it seems! Let me know if anything is unclear, if I made a glaring typo, or if you need further help.