Say $C$ is polyhedral cone, i.e., $C$ is written as a finite intersection of half spaces, $C=\{x \in \mathbb{R^n}: Ax\geq 0\}$. However, we can write $C$ also has a positive combinations of extreme rays, i.e., $C =\{Ry: y\geq 0\}$, where $R = [r_1,\dots, r_m] \in \mathbb{R}^{n\times m}$ is the juxtaposition of the extreme rays of $C$.
My question is then fairly simple: given matrix $A$, how do I compute matrix $R$?
The pair of matrices $(A, R)$ is called a Double Description Pair or DD Pair or Minkowski-Weyl Pair.
There is an algorithmic procedure known as the Double Description Method which given $A$, constructs $R$. The procedure is very similar to Fourier Elimination which is the oldest (and possibly the most inefficient) method for solving systems of linear inequalities.
This method roughly squares the number of columns of the computed $R_i$ in the $i$-th step resulting in exponentially big matrices $R$. Of course, not all of them are extreme rays and you can drop all but the extreme at each step but it's again possible that you will end up with an exponentially large number of extreme rays (compared to the size of $A$).
For more information check this paper and its references.
You may also want to check an implementation of this algorithm in Sage.
