I'm assuming you're thinking of some specific matrix representation $X_i \in \mathfrak{g}$ (let's assume it's the defining representation). Compute the Killing form, $\kappa_{ij} \doteq Tr (X_i\cdot X_j)$ (actually usually this is defined in the adjoint representation, but any faithful rep will do). The quadratic Casimir is then simply $ X_i \kappa_{ij} X_j$ (Einstein convention).
Other Casimirs can be obtained from the characteristic (secular) equation: Define $X(\omega) = \omega^i X_i$. The characteristic equation is
$$\det\left( X(\omega) - \lambda I \right) = \sum\limits_{j} (-\lambda)^{N-j} \phi_j(\omega) \equiv 0. $$
($N$ is the matrix dimension, and/or the dimension of the Lie algebra if you're using the adjoint representation). If you now perform the substitution $\omega^i \mapsto X_i$ in the coefficients $\phi_j (\omega)$, you get Casimir invariants $\phi_j (X)$!
It might seem like the higher the representation, the more invariants you get, but in fact all the invariants can be expressed in terms of the fundamental invariants of the defining representation. I can't think of many references right at the moment, but e.g. Gilmore: Lie groups, physics and geometry pp. 140 has a nice explanation. Also, google the boldface texts above.
