It's pretty straightforward if you use index notation.
A Brief Explanation of Index Notation
The rules for index notation are:
- summation is always implied when the same index occurs exactly twice in a single term. Such indices are called summing indices or dummy indices. Example: $\vec v \cdot \vec w=\sum_{i=1}^3v_i w_i = v_iw_i$. Here $i$ is the only summing index.
- if there is an index which is unpaired, then it is not summed over and must appear in every term on BOTH sides of the equation. Such indices are called free indices. Example: $x_a = R_{ab} v_b = \sum_{b=1}^3 R_{ab} v_b$. Here $a$ is the free index and $b$ is the summing index. Notice here that $x_a$ is the $a$th component of the vector $\vec x$ and $R_{ab}$ is a scalar indexed by $2$ indices -- where, in general, $R_{1,2} \ne R_{2,1}$ (this is how you denote the entries of a matrix in index notation).
- if your sum includes more than $2$ of the same index then you CANNOT use this implied summation convention. Example: $A_{ab}B_{cb}x_bv_c$ is not a valid expression using the implied summation convention because there are $3$ $b$ indices in this single term and thus the summation is ambiguous. If this summand occurs, you'd have to explicitly use the $\sum$ symbol.
A couple of symbols also need to be explained:
- The permutation symbol: $\epsilon_{pqr} = \begin{cases} 1, & pqr \in \{123, 231, 312\} \\ -1, & pqr \in \{321, 213, 132\} \\ 0, & \text{otherwise}\end{cases}$
- The Kronecker delta: $\delta_{pq} = \begin{cases} 1, & p=q \\ 0, & \text{otherwise}\end{cases}$
- The partial derivative wrt the $p$th argument of the function $f(x_1, x_2, \dots)$ will be denoted $\partial_p f$. For example, we'll denote $\dfrac{\partial f(x,y,z)}{\partial y}$ as $\partial_2 f$.
So now, let's start:
$$\begin{array}\ [\nabla \times (\nabla \times \vec F)]_p & \text{Just consider the $p$th component}
\\ = \epsilon_{pqr}\partial_q(\nabla \times \vec F)_r & \text{$[a\times b]_p := \epsilon_{pqr}a_qb_r$ -- this can be used as}
\\ & \text{the definition of the cross product}
\\ = \epsilon_{pqr}\partial_q\epsilon_{rst}\partial_sF_t & \text{Applying the same definition}
\\ = \epsilon_{pqr}\epsilon_{rst}\partial_q\partial_sF_t & \text{$\epsilon_{rst}$ is just a scalar and so it commutes}
\\ & \text{with derivative operator}
\\ = \epsilon_{rpq}\epsilon_{rst}\partial_q\partial_sF_t & \text{$\epsilon_{pqr}=\epsilon_{qrp}=\epsilon_{rpq}$}
\\ = (\delta_{ps}\delta_{qt} - \delta_{pt}\delta_{qs})(\partial_q\partial_sF_t) & \text{Apply the identity (which you can prove}
\\ & \text{if you want) $\epsilon_{rpq}\epsilon_{rst} = \delta_{ps}\delta_{qt} - \delta_{pt}\delta_{qs}$}
\\ = \partial_t\partial_pF_t -\partial_s\partial_sF_p & \text{Apply the identity $\delta_{qp}a_p = \delta_{pq}a_p=a_q$}
\\ = \partial_p\partial_tF_t - \partial_s\partial_sF_p & \text{Clairaut's theorem}
\\ = [\nabla (\nabla \cdot \vec F)]_p -[(\nabla \cdot \nabla)\vec F]_p & \text{We're almost done now}
\\ = [\nabla(\nabla\cdot \vec F)-\nabla^2\vec F]_p\ \ \ \ \ \ \square\end{array}$$
Once you get used to index notation, you'll see that there wasn't really any thinking required for this -- it was all pretty mechanical.
BTW, this process also shows you that $\vec F$ must have continuous second partials or else this identity doesn't hold (because we had to make use of Clairaut's theorem).