First observation: any combination of linear equality and inequality constraints can be converted to all $\le $ constraints. If you start with $$\begin{array}{c}
Ax\le a\\
Bx\ge b\\
Ex=e
\end{array},
$$you can convert it to the equivalent system of inequalities
$$
\left[\begin{array}{r}
A\\
-B\\
E\\
-E
\end{array}\right]x\le\left[\begin{array}{r}
a\\
-b\\
e\\
-e
\end{array}\right].
$$
Second observation (which requires knowledge of linear programming): if your feasible set is unbounded, you can slap an objective function on it that will improve along the ray (say, maximize the sum of the $x_i$, assuming $x\ge 0$) and get an unbounded linear program. Most LP solvers can find a ray once they have established that an LP is unbounded. If you prefer, you can try to apply the primal simplex method by hand. At some point you will encounter a basis where a variable wants to enter the basis (to improve the objective function) but there is no row in which to pivot. From that basic feasible solution you can easily identify a ray.
Addendum: How to use the all-negative simplex column. Assume a problem in "standard form" (meaning slacks/surpluses have been added) with constraints $Ax=b$. Let $B$ be the basis matrix when the negative column occurs, and for convenience assume the basic columns are the first few, so that $A$ gets partitioned as $[B\ \vert\ N]$ and your tableau corresponds to the system $I x_B + B^{-1}N x_N = B^{-1}b$. Let $k$ be the index of the all negative column and let $h$ be the column. If we freeze all the other nonbasic variables at 0, the system reduces to $x_B + h\cdot x_k = B^{-1}b$, with the right-hand side nonnegative. Since $h$ is all nonnegative, you can set $x_k$ as high as you want and the adjusted basic solution $x_B = B^{-1}b - h\cdot x_k$ will be nonnegative. So the ray originates at the current point $x_B = B^{-1}b, x_N = 0$ and moves in the direction $$\left[\begin{array}{c}
-h\\
1\\
0
\end{array}\right],$$ where I am assuming that $x_k$ is the first variable after the basic variables and the 0 has dimension equal to the number of nonbasic variables minus 1.
