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 $\pi:G\to G/H$ be the canonical coset projection. If the natural left action of $G$ on $G/H$ is not effective, replace $G$ by $G/N$, where $N\subset G$ is the closed, normal subgroup $N$ consisting of the elements that act trivially on $G/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 exists a unique non-degenerate, symmetric $n$-by-$n$ matrix of constants $Q$ such that $\pi^*g = Q_{ij}\,\omega^i\circ\omega^j$.
Since $Q$ is constant, there exists a unique $n$-by-$n$ matrix $\theta = (\theta^i_j)$, whose entries are left-invariant $1$-forms on $G$, such that
$$
\mathrm{d}\omega = {}-\theta\wedge\omega
\qquad\text{and}\qquad
Q\theta + {}^t(Q\theta) = \mathrm{d}Q = 0.
$$
(This is just the Fundamental Lemma of (pseudo-)Riemannian Geometry in this
context.) Note that applying $\mathrm{d}$ to both sides of this equation gives
$0 = \Theta\wedge\omega$, where, of course, $\Theta = \mathrm{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
$$
\mathrm{d}a = {}-\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 + (\mathrm{d}b +\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 = r^i_{jkl}\,a^k\omega^l$, where the $r^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
$$
\mathrm{d}b = -\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 $\mathrm{d}a$ and $\mathrm{d}b$ 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 $(G/H,g)$ is a (pseudo-)Riemannian space form, so, usually, one must differentiate these equations. The derivative of the $\mathrm{d}a$-equation won't give anything new, so one must differentiate the $\mathrm{d}b$-equation. This yields equations of the form
$$
0 = \mathrm{d}(\mathrm{d}b) = B,
$$
where $B = (B^i_j)$ and $B^i_j = (u^i_{jklm}a^m + v^{ip}_{jklq}b^q_p)\,\omega^k\wedge\omega^l$ for some explicit constants
$u^i_{jklm}= -u^i_{jlkm}$ and $v^{ip}_{jklq} = - v^{ip}_{jlkq}$.
It follows that the $a^i$ and $b^i_j$ are subject to the constant-coefficient linear relations
$$
u^i_{jklm}a^m + v^{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 $\mathrm{d}a$ and $\mathrm{d}b$ 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 that 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 $\mathrm{d}a_1$ and $\mathrm{d}a_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.
