I would like to study the characteristics of the non-linear dynamical system detailed below; in particular, I would like to find its set of fixed points; and to compute (i) the maximum values of particular variables (indicated below in red and blue) in the system as it evolves, and (ii) the time points at which these maxima are attained. I am looking for resources and advice to begin tackling these questions.
The dynamical system consists of a hierarchy of $K$ levels, where the $k$'th level consists of a vector $\mathbf{x}_k \in \mathbb{R}^{d_k}$ which is updated at each discrete time step $t$ as follows:
$$\mathbf{x}^{t}_k = \mathbf{x}^{t-1}_k \odot (W_k \cdot \left(\color{red}{ \mathbf{x}^{t}_{k-1} \oslash (\hat{W}_k^{\text{T}} \cdot \mathbf{x}^{t-1}_k)}\right) + \left( \color{blue}{ (\hat{W}_{k+1}^{\text{T}} \cdot \mathbf{x}^{t-2}_{k+1}) \oslash \mathbf{x}^{t-1}_k}\right))$$
The operations $\odot$ and $\oslash$ represent element-wise multiplication and division respectively. Because each vector $\mathbf{x}_k$ at a given level of the hierarchy is updated depending on the value of the vectors at other levels of the hierarchy $\mathbf{x}_{k-1}^t$ & $\mathbf{x}_{k+1}^{t-2}$, which might be of different dimensionality, each level has a matrix of constants $W_k \in \mathbb R^{d_k \times d_{k-1}}$ associated with it to transform a lower level vector to the current space or vice versa $\dagger$. The updates proceed in a particular temporal sequence, such that vectors at lower levels (e.g., level $k-1$) are updated prior to vectors at higher levels (e.g., level $k$).
There is an extra wrinkle: if any element in $\mathbf{x}_k^{t-1}$ is smaller than a pre-specified constant $\epsilon_1$, we replace that element with $\epsilon_1$ before computing $\mathbf{x}_k^{t}$; similarly, if any element in the dividend vector (in the elementwise divisions) is lower than $\epsilon_2$, we replace that element with $\epsilon_2$. So this is the actual update rule:
$$\small \mathbf{x}^{t}_k = \max(\epsilon_1,\mathbf{x}^{t-1}_k) \odot \left[(W_k \cdot \left(\mathbf{x}^{t}_{k-1} \oslash (\max(\epsilon_2,\hat{W}_k^{\text{T}} \cdot \mathbf{x}^{t-1}_k)\right) + \left( (\hat{W}_{k+1}^{\text{T}} \cdot \mathbf{x}^{t-2}_{k+1}) \oslash \max(\epsilon_2,\mathbf{x}^{t-1}_k)\right)\right] $$
From simulation "experiments" with one configuration of epsilons and $W$ matrices, I found that the multiplicative factor either approaches $0$ or $1$, which results in a stable system. However, this stability strongly depends on the choice of epsilon values. One reason I am interested in understanding this system is to have a principled method of finding epsilon values that keep the system stable, without having to resort to "grid search" simulation experiments.
$\dagger$ Note that $\hat W$ indicates a "normalized" $W$, where each row is divided by the sum of its elements.