I'm looking for attracting fixed points of the following differential equation in a vector $F$ of length $n$, and $M$ is a known $n \times n$ square matrix: \begin{align} \frac{dF}{dt} = I(F), \quad I(F) = (1-F) * (M \cdot F) - F * (M^{T} \cdot (1 - F)) \end{align} where $*$ means element-wise multiplication between vectors, and $\cdot$ means standard matrix-vector multiplication.
There is an additional set of constraints on $F$, since it physically represents a fraction of maximum occupancy. So each entry in the vector $F$ will be in the unit interval $[0, 1]$. Also, the dynamics will move things around but never create/destroy things, so $\sum_{i} F_{i} = n\tilde{c}$ is fixed. (This is manifest in the differential equation)
Time evolution with something like Runge-Kutta works in principle but is quite slow. The other thing I thought to try is a zero solver, but I'm having trouble with including the "fixed total occupancy" constraint in a way that doesn't make this equation way higher order than it needs to be. For example, adding a regularization term on the function \begin{align} L(F) = I(F) \cdot I(F) + \lambda \left( F \cdot 1 - n \tilde{c} \right)^{2} \end{align} which is zero iff both $I(F)$ is zero and the occupancy constraint is satisfied. However, this equation is now quartic in $F$.
Another trick I know is to let the last entry of $F$ be a function of the first $n-1$ entries, which would enforce the total occupancy constraint but would not address the problem of $F_{i} \in [0, 1]$.
Right now, I'm really not sure what the best approach would be to numerically approximate the fixed point, so any help would be greatly appreciated!