Probably, the easiest method is this (at least in characteristic zero, which I will assume henceforth): Suppose that $V\subset \mathbb{A}^n_k$ is the set of zeros of a polynomial ideal $I\subset k[x^1,\ldots,x^n]$, say, generated by some finite set $\{f_1,\ldots,f_m\}\subset k[x^1,\ldots,x^n]$. Let $\frak{a}$ be the (finite dimensional) Lie algebra of affine vector fields on $\mathbb{A}^n_k$ (i.e., the set of derivations of $k[x^1,\ldots,x^n]$ that preserve the subspace of polynomials of degree at most $1$). Consider the linear map $\Phi:{\frak{a}}\to k^m\otimes k[x^1,\ldots,x^n]/I$ defined as
$$
\Phi(X) = \bigl([X(f_1)]_I,\ldots,[X(f_m)]_I\bigr).
$$
The kernel of $\Phi$, say ${\frak{g}}_I\subset{\frak{a}}$, is the space of affine symmetry vector fields of the ideal $I$, and it is the Lie algebra of the connected (in the appropriate sense) subgroup $G^0_I$ of the set of affine transformations of $\mathbb{A}^n_k$ that preserve $I$, i.e., it is the irreducible component containing the identity of the algebraic group $G_I$ that consists of all affine symmetries of $I$. In particular, $G^0_I$ is an algebraic subgroup of the group of affine transformations of $\mathbb{A}^n_k$, and ${\frak{g}}_I$ is its (Zariski) tangent space at the identity. (Computing the full symmetry group $G_I$ is not a linear problem and is much harder in general, though one does know that $G^0_I$ is a normal subgroup of $G_I$ and that the quotient group $G_I/G^0_I$ is finite.)
In practice, if $I$ is complicated, computing $[X(p)]_I$ in an efficient way requires using Gröbner bases or some such tool. In all cases, though, the computation of ${\frak{g}}_I$ reduces to the computation of the kernel of a linear map between finite-dimensional vector spaces.
When $I$ is generated by a single polynomial $f$, one is just asking whether the 'remainder' of $X(f)$ 'divided by' $f$ is zero. (Here, I am intending the use of a so-called 'multivariate division algorithm' à la Buchberger's algorithm using a Gröbner basis, though a more naïve degree approach will work, too.) This is fairly easy to compute explicitly using a Gröbner basis, one, say, that uses some total order on the monomials in $k[x^1,\ldots,x^n]$, so the 'continuous' symmetries of a polynomial are not hard to compute using just about any symbolic algebra program.
For example, given $f$ of degree $d>0$, you can always choose coordinates $y^1,\ldots,y^n$ so that
$$
f = (y^n)^d + c_1(y^1,\ldots,y^{n-1})(y^n)^{d-1} + \cdots + c_d(y^1,\ldots,y^{n-1})
$$
where each $c_i$ is a polynomial of degree at most $i$. Clearly, there is a unique linear map $q:{\frak{a}}\to k$ such that
$$
X(f) = q(X)\,f + R(X,f)
$$
where $R(X,f)$ has degree at most $d$ in the $y^i$ and at most degree $d{-}1$ in $y^n$.
Now, $X(f)$ is a multiple of $f$ if and only if $R(X,f)=0$, and, for a given $f$, this is a finite number of linear equations for $X$.
