This sort of calculation is appropriately called index gymnastics. Here is the gory detail for calculating the covariant derivative of the vector field $\bar F$,
$$\begin{array}{rcll}
\displaystyle\frac{\partial \bar F}{\partial u^j}
&=& \displaystyle\frac{\partial (f^i\bar a_i)}{\partial u^j}
& \textrm{(expand in coordinate basis)} \\
&=& \displaystyle\frac{\partial f^i}{\partial u^j} \bar a_i
+ f^i \frac{\partial \bar a_i}{\partial u^j}
& \textrm{(chain-rule)} \\
&=& \displaystyle\frac{\partial f^i}{\partial u^j} \bar a_i
+ f^i \Gamma^k_{ij} \bar a_k
& \textrm{(definition of }\Gamma)\\
&=& \displaystyle\frac{\partial f^i}{\partial u^j} \bar a_i
+ f^k \Gamma^i_{kj} \bar a_i
& \textrm{(switch dummy indices)}\\
&=& \displaystyle\frac{\partial f^i}{\partial u^j} \bar a_i
+ f^k \Gamma^i_{jk} \bar a_i
& \textrm{(torsion-free)}\\
&=& \displaystyle\left(\frac{\partial f^i}{\partial u^j} \bar a_i
+ f^k \Gamma^i_{jk}\right) \bar a_i
& \textrm{(factor out }\bar a_i).
\end{array}$$
We have used your implicit assumption that the connection be torsion-free, $\Gamma^i_{jk} = \Gamma^i_{kj}$.
Of course, when you get used to doing such calculations you need not write so many steps.
Note that $i$ and $j$ in $a_{ij}b^{ij} = \sum_{i=1}^n \sum_{j=1}^n a_{ij}b^{ij}$ are dummy indices, they are being summed over.
(This is Einstein summation notation, a discovery Einstein was quite proud of, and rightfully so.)
We can change dummy indices so long as the index we change to is free.
For example, $a_{ij}b^{ij} = a_{ik}b^{ik}$, but $a_{ij}b^{ij} \ne a_{ii}b^{ii}$, since $a_{ii}b^{ii} = \sum_{i=1}^d a_{ii}b^{ii}$.