# Steady states and eigenvalues for a non-linear system:

I am trying to understand a non-linear system of equations, and find their steady states and dynamics. I am noob to understanding Mathematica (I am using version 6, but I have access to the latest versions in the lab), and I wondered if some kind soul would be able to help me with some of the code.

The system is below:

Clear[PBC];
PBC = {
d'[t] == (1/z) (α - β (d[t]/Y[t]) - γ L[t] - δ ψ[t]  - τ Y[t]) + ε (d[t]/ Y[t]) d[t],
L'[t] == (j (Y[t]/f) - e)/N,
ψ'[t] == p L[t] - \[Xi],
Y[t] == (v + b + α - β ( d[t]/Y[t]) - γ L[t] - δ ψ[t])
(ψ[t] + (1 - s) (1 - ψ[t]) - m z)/(1 - (1 - τ))
}
Clear[PBCSS];
PBCSS = Solve[PBC[[1 ;; 3]] /. {g_'[t] -> 0, g_[t] -> g}, {d, L, ψ}] // Simplify
Clear[PBCEigenSystem];
PBCEigenSystem = D[Transpose[PBCSS][[2, 1 ;; 3]], {{d, L, ψ}}] // Simplify) // MatrixForm


I know that I have done some no-no's (already edited in the code above I think), but the code seems to work. However, I am getting some problems with the output. For one, it does not solve all the system. I included some parameters for a numerical solution, but even if you take away the parametrization, it gives the same error as it does not solve for all variables (I have more choppy code to get some solutions for a large n number of iterations, I do not get it for steady states). Furthermore, I seem to be getting an error when I try to get the eigenvalues of the matrix. Also I have been trying to put a Manipulate handle on the steady states and eigenvalue estimation to see where the bifurcations are but I get errors.

Any help will be most appreciated!

-
Note that I is the built-in symbol for the imaginary number, which evaluates to a complex number: MatchQ[I + 1, _Real] returns False. Furthermore, as you already know it, E and even K have built-in meanings which might mess up your calculation. –  István Zachar Mar 22 '13 at 15:51
Thanks! Point taken, but even if all those parameters are replaced I still get the same errors. –  vyperultra Mar 22 '13 at 16:21
Are you sure your 4th equation is ok? Shouldn't it be the derivative, Y'[t] == ... instead of Y'[t] = ...? At the moment it is recursive. Stell there is no solution for steady state (Solve returns {}). –  István Zachar Mar 22 '13 at 17:20
Please give a minimal example, not your actual code. –  acl Mar 22 '13 at 17:40
@acl What's wrong with this code being not minimal (apart of possible errors)? –  István Zachar Mar 22 '13 at 17:50

Solve[PBC[[1 ;; 3]] /. {g_'[t] -> 0,
g_[t] -> g, α -> 0.05, γ -> 0.75, δ -> 0.75,
X -> 1, β -> 0.05, j -> 0.45, τ -> 0.35, v -> 1,
s -> 0.85, p -> 0.75, b -> 1, m -> 1.2, e -> 0.035, z -> 5.5,
f -> 3.5, N -> 4, ϵ -> 0.75, ξ -> 0.35}]


returns

{{ψ -> -0.527037 - 0.244898 d + 20.2041 d^2, L -> 0.466667, Y -> 0.272222}}


so {d, L, ψ} are not suitable variables to solve for. In fact, L and Y turn out to be constants and d and ψ are not independent.
I think you have to rework the math of your problem.

-