3027.sagews

T, IT, V_I, V_NI, s, p, Tmax, Dt, k, delta, N, c, n_RT, n_p = var('T IT V_I V_NI s p Tmax Dt k delta N c n_RT n_p')

dT = s + p*T*(1-(T+IT)/Tmax) - Dt*T- (1-n_RT)*k*V_I*T
dT

\newcommand{\Bold}[1]{\mathbf{#1}}{\left(n_{\mbox{RT}} - 1\right)} T V_{I} k - {\left(\frac{\mbox{IT} + T}{\mbox{Tmax}} - 1\right)} T p - \mbox{Dt} T + s

dIT = (1-n_RT)*k*V_I*T - delta*IT
dIT

\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(n_{\mbox{RT}} - 1\right)} T V_{I} k - \mbox{IT} \delta

dV_I = (1-n_p)*N*delta*IT - c*V_I
dV_NI= n_p*N*delta*IT - c*V_NI
dV_NI

\newcommand{\Bold}[1]{\mathbf{#1}}\mbox{IT} N \delta n_{p} - V_{\mbox{NI}} c

J = matrix([[diff(dT,T), diff(dT,IT), diff(dT,V_I), diff(dT,V_NI)],[diff(dIT,T), diff(dIT,IT), diff(dIT,V_I), diff(dIT,V_NI)], [diff(dV_I,T), diff(dV_I,IT), diff(dV_I,V_I), diff(dV_I,V_NI)], [diff(dV_NI,T), diff(dV_NI,IT), diff(dV_NI,V_I), diff(dV_NI,V_NI)]])

J1 = J.subs(fixedpoints[0])
J2= J1.subs(delta=0.24);
J3= J2.subs(k=2.4e-5);
J4=J3.subs(c=2.424);
J5= J4.subs(T=500)
J6=J5.subs(N=202)
J7=J6.subs(Tmax=1500)
lambdas1 = J7.eigenvalues()
lambdas1

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{3}{250} \, \sqrt{4040 \, {\left(n_{p} - 1\right)} n_{\mbox{RT}} - 4040 \, n_{p} + 12321} - \frac{333}{250}, \frac{3}{250} \, \sqrt{4040 \, {\left(n_{p} - 1\right)} n_{\mbox{RT}} - 4040 \, n_{p} + 12321} - \frac{333}{250}, -\frac{303}{125}, -\mbox{Dt} + \frac{1}{3} \, p\right]

J1 = matrix([[diff(dIT,IT), diff(dIT,V_I), diff(dIT,V_NI)], [diff(dV_I,IT), diff(dV_I,V_I), diff(dV_I,V_NI)], [diff(dV_NI,IT), diff(dV_NI,V_I), diff(dV_NI,V_NI)]])
lambdas= J1.eigenvalues()
lambdas

\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, c - \frac{1}{2} \, \delta - \frac{1}{2} \, \sqrt{4 \, {\left({\left(\delta k n_{p} - \delta k\right)} N n_{\mbox{RT}} - {\left(\delta k n_{p} - \delta k\right)} N\right)} T + c^{2} - 2 \, c \delta + \delta^{2}}, -\frac{1}{2} \, c - \frac{1}{2} \, \delta + \frac{1}{2} \, \sqrt{4 \, {\left({\left(\delta k n_{p} - \delta k\right)} N n_{\mbox{RT}} - {\left(\delta k n_{p} - \delta k\right)} N\right)} T + c^{2} - 2 \, c \delta + \delta^{2}}, -c\right]

fixedpoints = solve([dIT == 0, dV_I == 0, dV_NI == 0], [IT, V_I,V_NI], solution_dict=True)
fixedpoints

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{V_{\mbox{NI}}:\: 0, V_{I}:\: 0, \mbox{IT}:\: 0\right\}\right]