Is there a specific formula/method to find geodesics for a Homogeneous space? (excluding general methods applicable to arbitrary riemannian manifold)

I assume you mean riemannian homogeneous space and you are talking about geodesics relative to the LeviCivita connection. If so, you can always try to find geodesics which are homogeneous; that is, geodesics which correspond to the orbits of oneparameter subgroups. Let $M = G/H$ and let $\mathfrak{g} = \mathfrak{h} \oplus \mathfrak{m}$ be a reductive split (always possible for riemannian homogeneous spaces, since $H$ is compact) of the Lie algebra of $G$. By homogeneity, you can look for geodesics passing through the identity coset $o \in M$ and then transport them using the $G$action to any other point. To specify a geodesic through $o$ it is necessary and sufficient to specify the tangent vector at $o$, which we can think of as element $X \in \mathfrak{g}$. Then the orbit $\exp(t X) \cdot o$ is a geodesic if and only if $$ \left< [X,Z]_{\mathfrak{m}},X_{\mathfrak{m}} \right> = 0 $$ for all $Z \in \mathfrak{m}$, and where the subscript means the component along $\mathfrak{m}$. Although this is not guaranteed to find you all geodesics, there are some homogeneous spaces (socalled geodesic orbit or g.o. spaces) for which all geodesics are of this type. 


One method is to use the generalization of the principle of conservation of momemntum and angular momentum. If $\phi_t$ is any one parameter group of isometries of any Riemannian manifold, with time derivative expressed as a vector field $X$, then for any geodesic $g$, the inner product $\left \langle X, \dot g \right \rangle $ is constant along the geodesic, where $\dot g$ denotes the unit tangent vector to $g$ (assuming $g$ is parametrized by arclength). In Euclidean space, applied to 1parameter groups of translations, this gives the principle of conservation of momemntum for a geodesic; applied to 1parameter groups of rotations, it gives the principle of conservation of angular momentum. In a homogeneous space, there are at least enough such vector fields to form a basis at each point. In an $n$dimensional homogeneous space, if there are more than $n$ linearly independent infinitesimal isometries, then the equations, for a typical initial vector, are consistent only in a proper subset, and on that subset, they typically define a vector field. This may give you all the information you need, without even integrating the vector field. But in any case, it cuts the differential equation down to many fewer degrees of freedom. 
