Here is an algorithm to compute the Lie algebra of the group of isometries of a homogeneous space $G/H$ endowed with a $G$-invariant (pseudo-)Riemannian metric $g$. It is phrased in terms of essentially algebraic computations using the left-invariant forms on $G$, but it could be reduced completely to computations with the Lie algebra (and the matrix $Q$ that defines the metric) if that is what one wanted to do.
Let $\frak{g}$ and $\frak{h}\subset\frak{g}$ denote the Lie algebras of $G$ and $H\subset G$ respectively. Set $s = \dim\frak{h}$ and let $n>0$ be the dimension of $\frak{m} = \frak{g}/\frak{h}$.
Let $\omega = (\omega^i)$ (where $1\le i\le n$) be a basis for the left-invariant $1$-forms on $G$ such that $\omega = 0$ defines the foliation of $G$ by the left cosets of $H$. Then there will exist a unique non-degenerate, symmetric $n$-by-$n$ matrix $Q$ such that the quadratic form $Q_{ij}\omega^i\circ\omega^j$ is the $\pi$-pullback to $G$ of $g$ under the natural coset mapping $\pi:G\to G/H$.
Since $g$ is invariant under the action of $G$, there exists a unique $n$-by-$n$ matrix $\theta = (\theta^i_j)$, whose entries are left-invariant $1$-forms on $G$, such that
$$
d\omega = {}-\theta\wedge\omega
\qquad\text{and}\qquad
Q\theta + {}^t(Q\theta) = 0.
$$
(This is just the Fundamental Lemma of (pseudo-)Riemannian Geometry in this
context.) Note that applying $d$ to both sides of this equation gives
$0 = \Theta\wedge\omega$, where, of course, $\Theta = d\theta + \theta\wedge\theta$, is the curvature of the connection $\theta$ and hence can be written in the form $\Theta = R(\omega\wedge\omega)$, where the coefficients in $R$ are constants, since $R$ is left-invariant as a function on $G$.
Suppose now that a vector field $Y$ on $G$ be $\pi$-related to a $g$-Killing vector field $Z$ on $G/H$, and let $\omega(Y) = a$. Then the $g$-Killing equation for $Z$ implies that
$$
da = {}-\theta\ a + b\ \omega
$$
where $b$ is an $n$-by-$n$ matrix of functions that satisfies $Qb + {}^t(Qb)=0$.
Taking the exterior derivative of this equation gives
$$
0 = {} -\Theta\ a + (db +\theta\ b - b\ \theta)\wedge\omega.
$$
Now, by counting dimensions to show that the corresponding linear algebra problem always has a unique solution, it is easy to see (and, in any given case, explicitly compute) that there exists a matrix $\rho = (\rho^i_j)$ of $1$-forms such that $\Theta\ a = \rho\wedge\omega$ and such that $Q\rho + {}^t(Q\rho) = 0$. In fact, one has $\rho^i_j = c^i_{jkl}a^k\omega^l$, where the $c^i_{jkl}$ are constants determined by the curvature form $\Theta$. (The exact formula is not important for the following argument.) Thus, the above equation can be written as
$$
db = -\theta\ b + b\ \theta + \rho(a,\omega),
$$
where I have written the term $\rho$ as $\rho(a,\omega)$ to emphasize that this is some constant-coefficient bilinear pairing of $a$ and $\omega$ taking values in the Lie algebra ${\frak{so}}(Q)$.
The above equations for $da$ and $db$ are then a total (linear) differential system whose solutions give the Lie algebra of $g$-Killing vector fields on $G/H$.
Now, this system is not Frobenius unless the $(G/H,g)$ is a (pseudo-)Riemannian space form, so you have to differentiate these equations. The derivative of the $da$-equation won't give anything new, so one must differentiate the $db$-equation. This yields equations of the form
$$
0 = d(db) = B,
$$
where $B = (B^i_j)$ and $B^i_j = (b^i_{jklm}a^m + c^{ip}_{jklq}b^q_p)\omega^k\wedge\omega^l$ for some explicit constants
$b^i_{jklm}= -b^i_{jlkm}$ and other constants $c^{ip}_{jklq} = - c^{ip}_{jlkq}$. (I can't see anything wrong with the TeX here, but it's not typesetting correctly.)
It follows that the $a^i$ and $b^i_j$ are subject to the constant-coefficient linear relations $b^i_{jklm}a^m + c^{ip}_{jklq}b^q_p = 0$ in addition to the linear relations on $b$ already known: $Qb + {}^t(Qb) = 0$. Differentiating these new linear relations and using the $da$ and $db$ formulae to express the results in terms of left-invariant forms with coeffcients that are constant linear combinations of the $a$- and $b$-components, one might get further constant coefficient linear relations among the $a$- and $b$-components. Repeat this process with the new relations (if any) until no new linear relations are found.
At this point, the linear relations between the $a$- and $b$-components will have some solution space of dimension $n{+}s{+}r$ for some $r\ge0$ (but, necessarily, $r\le n(n{-}1)/2 - s$). It will then follow that the space of $g$-Killing vector fields on $G/H$ has dimension $n{+}s{+}r$. Moreover, one can compute the Lie algebra structure on this space by using the formula
$$
\omega\bigl([Y_1,Y_2]\bigr)
= Y_1(a_2) - Y_2(a_1) - \theta\wedge\omega\bigl(Y_1,Y_2\bigr)
$$
and the formulae for $da_1$ and $da_2$. Thus, the algebra structure of the $g$-Killing fields will follow directly by algebraic operations from the structure of the algebras $\frak{g}$ and $\frak{h}$ and $Q$.
In any given instance, this algorithm can be implemented on a computer without difficulty.
