Basically you are missing a term. Think of Lagrangian mechanics. If $q(t)$ is your generalized coordinate, and you make the change $q(t) \to q(t) + \delta q(t)$, then remember that the change in the action is $\delta S = \int_{t_i}^{t_f} \left(\partial_q L - \partial_t \partial_\dot{q} L \right) \delta q$. Thus you can basically think of the change in the Lagrangian $\frac{\delta L}{\delta q}$ to be $\partial_q L - \partial_t \partial_\dot{q} L$.
An analogous thing happens when you try to compute $\frac{\delta}{\delta\mathbf{A}}|\mathbf{B}|^2$. You can think of this derivative as $\partial_{A_i} B_j B^j - \partial_k \partial_{\partial_k A_i} B_j B^j$. (You should try to derive this expression for yourself. The steps are analogous to what you have to do to derive the euler-lagrange equations. In both, the key step is integration by parts.) In your case, the first term is zero, but the second term is not. Plugging in the definition of $\mathbf{B}$, we find that the second term is (including the minus sign) $-2 \partial_k \epsilon_{j}{}^{ki}\epsilon^{j}{}_{hl}\partial^h A^l = 2 \epsilon^{ik}{}_{j} \partial_k \epsilon^{j}{}_{hl} \partial^h A^l $. In vector form, this last expression is $2 \mathbf{\nabla} \times \left( \mathbf{\nabla} \times \mathbf{A} \right)$, which is hopefully what you are supposed to get.
Appendix: Why is $\frac{\delta}{\delta\mathbf{A}}|\mathbf{B}|^2 = \partial_{A_i} B_j B^j - \partial_k \partial_{\partial_k A_i} B_j B^j$
Review of derivation of Euler-Lagrange equations
Remember that in lagrangian mechanics, to get the equation of motion for a system from the lagrangian $L(q,\dot{q})$, you do the following. First you define the action $S$ for a path $q(t)$ to be $S[q]=\int L(q,\dot{q})\,dt$. Next you say that a physically allowed path is one that extremizes the action. That is, if $q_p$ is a physically allowed path, then $\frac{\delta S}{\delta q} \Big|_{q=q_p} =0$. Another way of saying this is that if you take the physical path $q_p$ and change it by any small amount $\delta q$, the action $S[q_p + \delta q]$ is the same as the action you started with $S[q_p]$ up to first order in $\delta q$.
To see what this actually means, plug the definition of $S$. You find that $$S[q_p + \delta q] = \int L(q_p + \delta q,\dot{q}_p + \delta\dot{q}) \, dt.$$ Taylor expanding $L$, you find $$S[q_p + \delta q] = \int L(q_p) + \partial_q L\, \delta q + \partial_\dot{q} L\, \delta \dot{q} \, dt.$$ The next step is the crucial step. You use integration by parts on the third term to move the time derivative from $\delta q$ to $\partial_\dot{q} L$. You then get $$S[q_p + \delta q] = \int L(q_p) + \partial_q L\, \delta q - d_t \partial_\dot{q} L\, \delta q\, dt.$$ Recognizing that the first term just gives you $S[q_p]$ and rearranging, you finally get $$S[q_p + \delta q]= S[q_p] + \int \left(\partial_q L-d_t \partial_\dot{q} L\right)\delta q \, dt.$$
What does this mean? Well we said that for any $\delta q$, $S[q_p + \delta q]$ should be the same as $S[q_p]$ up to first order in $\delta q$. This means, looking at the above equation, that $\int \left(\partial_q L-d_t \partial_\dot{q} L\right)\delta q \, dt=0$, no matter what $\delta q$ is. The only way this is possible is if $\partial_q L-d_t \partial_\dot{q} L=0$ for all $t$. This is how the Euler-Lagrange equation is found in lagrangian mechanics.
Two small modifications
Hopefully the derivation above wasn't new to you. But already it is essentially all you need for your problem. There are just two small modifications.
The first modification is if $q$ were a vector (or tensor or spinor). For example, $q$ could be a position in three dimensional space. Our derivation remains exactly the same except now $\partial_q$ is a gradient operator (let's call it $\nabla_q$, and similarly so is $\partial_\dot{q}$ (let's call it $\nabla_\dot{q}$). Then the equation $\partial_q L-d_t \partial_\dot{q} L=0$ becomes the vector equation $\nabla_q L-d_t \nabla_\dot{q} L=\mathbf{0}$, with the $\mathbf{0}$ on the right hand side being the zero vector.
Now let's go back and do another modification. Let's suppose the argument of $q$, which had been $t$, is now a vector instead of just a number. An example would be if $q$ represented the height of a two dimensional elastic sheet: the height not only depends on time (since the sheet is moving), but also depends on where on the sheet you are. Since $t$ is now a multidimensional object the derivative of $q$ with respect to $t$ (i.e., $\dot{q}$) is now a gradient and should be written as $\nabla_t q$. The lagrangian is then written as $L(q,\nabla_t q)$. Since $\nabla_t q$ is itself a vector, the partial derivative $\partial_\dot{q}$ becomes a gradient operator $\nabla_{\nabla_t q}$. The term $\partial_\dot{q} L \, \delta \dot{q}$ then becomes $\nabla_{\nabla_t q} L \cdot \delta \nabla_t q$. After the integration by parts step, this terms becomes $-\nabla_t \cdot \nabla_{\nabla_t q} L \, \delta q$. So the Euler-Lagrange equation becomes $\partial_q L - \nabla_t \cdot \nabla_{\nabla_t q} L = 0$.
Now you can also make both of these modifications at once. In fact that is the situation of your question. But at this point things get messy enough that you have to use index notation. I will first rewrite the three version of the Euler-Lagrange equation using index notation, then I will write the fourth, combined version. The first version doesn't have any indices. It is just $$\partial_q L - \partial_t \partial_\dot{q} L=0.$$
For the second version, we assume that $q$ is a multidimensional object with components $q_\alpha$. (Here I am using greek letters to index the components because I will use latin letters to index the field argument components, and these indices may no be the same. For example the field could be a spinor-valued.) In this case the Euler_Lagrange equation is easy to write, you just add in the $\alpha$ index: $$\partial_{q_\alpha} L - \partial_t \partial_{\dot{q}_{\alpha}} L=0.$$
Now for the third version, were $t$ is a vector with components $t_k$. As a shorthand, I will write the derivative with respect to $t_k$ simply as $\partial_k$ (instead of $\partial_{t_k}$). Then $\dot{q}$ becomes $\partial_k q$, and the Euler-Lagrange equation becomes $$\partial_q L - \partial_k \,\partial_{\partial_k q} L=0.$$
Now combining these two, we get the equation $$\partial_{q_\alpha} L - \partial_k \partial_{\partial_k\, q_\alpha} L=0.$$
Now it so happens that in your example your $\alpha$ index is not a spinor index but a plain old vector index, so it is the same type as your $k$ index, which can be confusing. But since it is the same type I used the latin index $i$. Also your lagrangian is $B_jB^j$, and your field "$q$" is the vector potential $A$. Plugging these in, we get the expression $$\partial_{A_i} B_jB^j - \partial_k \partial_{\partial_k\, A_i} B_jB^j=0.$$
