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}$.