Dear @tomleslie, @Carl Love, et al.

May I ask a conceptual follow-up question related to the numerical exercise above?

Consider the following two examples: Example1.mw and Example2.mw. Both examples are generated from the same linear system (of 4 equations in 4 unknowns) of the following type:

(EQ1): y[1]=f(y[1], y[2], y[3], y[4])

(EQ2): y[2]=f(y[1], y[2], y[3], y[4])

(EQ3): y[3]=f(y[1], y[2], y[3], y[4])

(EQ4): y[4]=f(y[1], y[2], y[3], y[4])

The difference between the two examples is that, for the parameter values in **Example1** the iterative process described above **converges**, while for parameter values in **Example2** the iterative process **diverges**.

Importantly, in both cases, there exists an (up-to-scale unique) positive real solution to the above system of equations (referred henceforth as an equilibrium), and in Example1 we can converge to it iteratively, while in Example2 we cannot.

In view of the above, my **question** is: Does the fact that the solution in Example2 cannot be reached iteratively reveal something about the **stability of the equilibrium** in this case? Simply put, is equilibrium in Example2 unstable and, if so, how do I test this in Maple?

I should also mention the following piece of information that may be relevant to help you think about my question: The above system of equations can be written in the following **matrix form**:

y=(I-A)^{-1} * x,

where y is a column vector of unknowns, I is an identity matrix, A is square matrix of independent variables and x is a column vector of independent variables (in my model, these are the "shocks" that I feed into the system to study their effect on the vector of y's for given values of A).

For the parameter values underlying **Example1**, the **spectral radius** (rho) of matrix A is smaller than one, i.e., **rho(A)<1**. In this case, matrix (I-A)^{-1} can be rewritten as Neumann series and *n*th element (degree) of this power series corresponds to the value of the unknowns after *n* iterations of the above system of linear equations. In contrast, for the parameter values underlying **Example2**, the spectral radius of A is larger than 1, **rho(A)>1**, which, obviously, precludes the Neumann series decomposition of this matrix. So, to paraphrase my question in matrix notation, does the case of **rho(A)>1 **reveal something about the stability of the equilibrium vector y and, if so, how do I test the **stability** in Maple?

Thank you very much in advance!