This question is related to this one Simple question on symmetric tensors.
To prove that a vector field $Z$ is Killing, we use the identity
$$0=(L_Zg)(X,Y)=g(X,\nabla_YZ)+g(\nabla_XZ,Y)\ \ \ \forall \ \ X,Y$$
It is clear that $(L_Zg)(X,Y)$ is a symmetric tensor. My Question is: Is it enough to prove that
$$(L_Zg)(X,X)=0\ \ \ \forall \ \ \ \ X$$
In this case we can say that $Z$ is a Killing vector field if
$$0=g(\nabla_XZ,X)\ \ \ \forall \ \ X$$
Note
The first identity has three solutions in the plane namely $(1,0),(0,1),(y,-x)$. But the second identity has a more general solution namely $(f(y),h(x))$. So I think they are not equivalent!!
First, let's compute the mentioned example in detail:
Example Find the Killing fields $Z$ of the standard metric $\bar{g} := dx^2 + dy^2$ on the plane $\mathbb{R}^2$.
If we write $$Z = f \partial_x + g \partial_y,$$ for some functions $f,g$ of $(x, y)$, then $$g(Z, \cdot) = f \,dx + g \,dy$$ and so the bilinear form $(X, Y) \mapsto g(\nabla_X Z, Y)$ is $$g(\nabla_{\cdot} Z, \cdot) = \nabla(g(Z, \cdot)) = (f_x \,dx + f_y \,dy) \otimes dx + (g_x \, dx + g_y \, dy) \otimes dy.$$ Now, inserting $X = p \partial_x + q \partial_y$, where $p, q$ are functions of $(x, y)$ into both slots gives
$$p^2 f_x + 2pq(f_y + g_x) + q^2 g_y.$$
Now, this must hold separately for every value of $p, q$, and so $f_x, f_y + g_x, g_y$ must all vanish separately. The first gives that $f(x, y) = j(y)$ for some constant $a$ and the third $g(x, y) = k(x)$. Substituting these in the remaining equation and rearranging gives $j_y = - k_x$. Now, the LHS depends only on $y$ and the RHS only on $x$, so both sides are actually constant, say, $c$. Substituting gives that the general solution is
$$Z = (c y + a)\partial_x + (-c x + b)\partial_y = a \partial_x + b \partial_y + c(y \partial_x - x \partial_y).$$
which is exactly the solution space claim.
As in my comment, I mentioned that it is insufficient to evaluate the quadratic expression on a basis to determine a full set of component conditions. If we substitute in $\partial_x$, we get the first condition $f_x = 0$ above and substituting $\partial_y$ gives $g_y = 0$. The solution to this pair of equations is indeed $(j(y), k(x))$ but we can see that this computation misses the critical cross term that cut us down to a finite-dimensional solution space.
Example Consider the symmetric bilinear form $Q = dx \,dy$ on a coordinate patch of a $2$-manifold. Then, $Q(\partial_x, \partial_x) = 0$ and $Q(\partial_y, \partial_y) = 0$, but $Q$ is not the zero symmetric bilinear form---so, it's insufficient to test a basis. Put another way, there is more than one quadratic function $q(x, y)$ on $\mathbb{R}^2$ whose graph contains the $x$- and $y$-axes.
As for the proof that checking the condition $g(\nabla_X Z, X) = 0$ is sufficient, your post essential includes it. A little more formally, one might say that if that condition holds for all $X$, then polarizing the symmetric $2$-tensor form $X \mapsto g(\nabla_X Z, X)$ gives (up to a nonzero constant factor) that the Killing equation holds.