The most difficult part is to actually get a set of consistent boundary conditions in the first place - this requires a combination of educated guessing, physical insights, prior experience with related problems, detailed calculations and trial-and-error. In short, it is a bit of an art.
However, once you have a set of boundary conditions (as in your case the NHEK boundary conditions) things are fairly straightforward.
Let me denote the asymptotic background metric by $g$ and the state-dependent fluctuations by $h$, so that any metric of the form $g+O(h)$ is allowed by the boundary conditions.
Your goal is to check what are the gauge transformations that preserve your boundary conditions. In pure gravity, these must be some diffeomorphisms generated by some vector field $\xi$, such that
\begin{equation}
L_\xi (g+h) = O(h)
\end{equation}
where $L_\xi$ is the Lie-derivative. Since this is an equation for a symmetric tensor in $D$ dimensions you get $D(D+1)/2$ independent linear first order PDEs for the vector field $\xi$.
If you tried to solve the equation above you were precisely doing the correct thing, which hopefully answers part of your question. Let me now address the other part, namely how to solve these PDEs.
In many examples you can solve for the most general vector field $\xi$ compatible with the condition above simply by guessing a suitable Ansatz and then showing that it works.
In most applications you have some series expansion of the asymptotic metric in powers of some radial coordinate $r$ (or in exponentials of $r$, depends on your gauge choice for the radial coordinate).
To reduce clutter let me assume that the various tensor components of $g$ and $h$ are expressed in some Laurent series of $r$, and that $r\to\infty$ corresponds to the asymptotic boundary.
Then you just make the same kind of power series Ansatz for the vector field $\xi$. Typically, all components of the vector field start at $O(1)$ or smaller, but this need not be the case. If you have no clue at all, then just make the Ansatz
\begin{equation}
\xi^0 = r^{n_0} (\xi^0_0 + \xi^0_1/r + ...)
\end{equation}
where the coefficient functions $\xi^0_i$ are allowed to depend on all boundary coordinates a priori. You make a similar Ansatz for all the other components of the vector field.
Evaluating the conditions from the Lie variation then determines the exponents $n_i$ and may put restrictions on the functions $\xi^i_j$.
See for instance exercise (17.1) in exercises of week 7 on my teaching webpage http://quark.itp.tuwien.ac.at/~grumil/teaching.shtml for a guided tour through the standard AdS$_3$ example. If you have never done this sort of calculation before I recommend you start with that example before breaking the NHEK.
[On a sidenote, I have never tried to implement this algorithm in Mathematica, since I usually work in 3 dimensions where a calculation by hand is fairly quick, but I do not see any reason why it should not work in Mathematica.]
