Consider first a one-dimensional dynamical system given by $dx/dt = f(x)$. Suppose that the origin is a fixed point, i.e. $f(0)=0$. Suppose that we're interested in whether trajectories that start from positive initial values of $x$ (hereafter "positive trajectories") will converge towards the origin or diverge away from it.
If $df/dx$ is negative at 0 then the fixed point is locally asymptotically stable; if it's positive then it's unstable. If $df/dx=0$ at the origin then the fixed point is called "non-hyperbolic". If $d^2 f/dx^2 <0$ then it may or may not be Lyapunov stable (because negative trajectories might diverge), but positive trajectories will converge toward the origin, albeit more slowly than a linearly stable fixed point. In general positive trajectories will converge towards the origin iff the lowest non-zero derivative is negative.
My question is about how to generalise this to the $n$-dimensional case. I'm dealing with systems where no trajectory starting in the non-negative orthant can leave the non-negative othant, and I'm interested only in trajectories that are confined to it. ("Nonnegative trajectories".)
In the case where the origin is a hyperbolic fixed point it's clear how to proceed. One creates the linear approximation $\mathbf{\dot{x}} = A\mathbf{x}$. The Jacobian matrix $A$ is the analog of $dx/dt$ - if all its eigenvalue have real part strictly less than 0 then the fixed point is stable, and hence all trajectories (positive or otherwise) will converge toward the origin. If one of the eigenvalues has real part greater than $0$ then the fixed point is unstable. (The condition that trajectories are confined to the non-negative orthant implies that $A$ can have negative entries only on its diagonal, which in turn means that its leading eigenvalue will always have a corresponding eigenvector in the nonnegative orthant, by the Perron-Frobenius theorem applied to $e^{\varepsilon A}$ for sufficiently small $\varepsilon$; hence we don't need to worry about the possibility that it might be unstable only for non-nonnegative trajectories.)
However, I'm currently dealing with systems where the leading eigenvalue of $J$ is exactly 0, i.e. the fixed point is non-hyperbolic. In this case I guess I have to move up to a quadratic approximation: $$ \frac{dx_i}{dt} = \sum_j A_{ij}x_j + \sum_{jk}B_{ijk}x_jx_k. $$ I can calculate the values of $A_{ij}$ and $B_{ijk}$ for my systems. (They are generated computationally, and in my systems $n$ is of the order 1000 or so.) In addition to $A$ having negative entries only on its diagonal, the condition of being confined to the positive orthant implies that if $B_{ijk}$ is non-zero then it is negative if and only if $i=j$ or $i=k$ (or both).
The problem is that I don't know how to tell from these values whether the system is stable or not. I'm looking both for an algorithm to determine this, and for necessary and sufficient conditions that are easy to compute. It strikes me that there's probably a nicely worked out theory that deals with this, but I don't know its name, so I can't find any information about it.