Given the angular momentum operators $$L_x,L_y,L_z, \mathbf{L}^2 \stackrel{\text{by def}}{=} \mathbf{L}\cdot \mathbf{L} = L_xL_x+L_yL_y+L_zL_z$$ satisfy the following commutation relations:$$[L_x,L_y] = i\hbar L_z,[L_x,L_z] = -i\hbar L_y, ... $$ which can be summarized using the Levi-Civita symbol: $$[L_i,L_j] = i\hbar\sum_{k}\epsilon_{ijk}L_k \stackrel{\text{Einstein convention}}{\longrightarrow} =i\hbar \epsilon_{ijk}L_k \tag{1}$$ Now it often proves to be useful to use the following commutation relation: $$ [L_i,\mathbf{L}^2] = 0, \forall i = 1,2,3$$ The standard way of proving this in Physics books (for example Sakurai or Townsend) is to just pick one of the components, for instance, $L_z$, and directly compute the commutator $$[L_z, L_x^2+L_y^2+L_z^2] = \cdots = 0$$ using $(1)$ and the well-known commutation identity for $[A,BC]$.
I wanted to compute this using only $(1)$ but I'm having trouble with the levi-civita symbols (I'm not at all proficient, unfortunately ): $$ \sum_{j=1}^{3}[L_i, L_j^2]= \sum_{j}[L_i,L_j]L_j+L_j[L_i,L_j] = \sum_{j}\left[i\hbar \left(\sum_k\epsilon_{ijk}L_k\right)L_j + i\hbar L_j\left( \sum_k\epsilon_{ijk}L_k \right) \right] = $$ $$ i\hbar \sum_{jk} \epsilon_{ijk}(L_kL_j+L_jL_k) = i\hbar \sum_{jk}\epsilon_{ijk}\left(2L_kL_j-[L_k,L_j] \right) = $$ $$ i\hbar\left( 2\sum_{jk}\epsilon_{ijk}L_kL_j-i\hbar \sum_{ijk}\epsilon_{ijk}\epsilon_{kji}L_i \right) = $$ $$\stackrel{{\epsilon_{kji} = -\epsilon_{ijk}}}{=} 2i\hbar\sum_{jk}\epsilon_{ijk}L_kL_j-\hbar^2\sum_{ijk}\epsilon_{ijk}\epsilon_{ijk}L_i $$ This is where I've been stuck: I can't think of a way that would lead to those two expressions ending up canceling each other out. Hence I'm here to see if anyone could either help me spot the mistakes I made or maybe unravel the final expression with the L-C symbols. Thanks in advance.
Look at $$\sum\limits_{j,k} \epsilon_{ijk} (L_j L_k+L_k L_j)\tag{1} \label{eq1}$$ in the second line of your calculation. $\epsilon_{ijk}$ is antisymmetric with respect to $j \leftrightarrow k$, whereas the anticommutator $\{L_j ,L_k\}:= L_j L_k +L_k L_j$ is symmetric under the same interchange of the indices $j,k$. As you sum over these indices in \eqref{eq1}, the resulting expression has to vanish.