x

y

flip2 = 0
ymin = -141.0889465212123
ymax = 171.2537462477026
xmin = -176.36308924798803
xmax = 177.97307176384336
dy = 10
dx = 10
M2 = 150
flip1 = 0
M1 = 200
theta2 = 100
theta1 = 100
sigma1 = 30
S0(x) = 1 / (1 + ℯ^((-x)))
sigma2 = -10
S1base(x) = (M1 * S0((x - theta1) / sigma1))
S2base(x) = (M2 * S0((x - theta2) / sigma2))
LL = Corner(1)
UR = Corner(3)
S1fs0 = CurveCartesian(s, theta1, s, xmin, xmax)
S1f = CurveCartesian(S1base(s), s, s, ymin, ymax)
S1nfs0 = CurveCartesian(theta1, s, s, ymin, ymax)
S1nf = CurveCartesian(s, S1base(s), s, xmin, xmax)
S1 = If(flip1 ≟ 1, If[sigma1 ≟ 0, S1fs0, S1f], If[sigma1 ≟ 0, S1nfs0, S1nf])
S2fs0 = CurveCartesian(s, theta2, s, xmin, xmax)
S2f = CurveCartesian(S2base(s), s, s, ymin, ymax)
S2nfs0 = CurveCartesian(theta2, s, s, ymin, ymax)
S2nf = CurveCartesian(s, S2base(s), s, xmin, xmax)
S2 = If(flip2 ≟ 1, If[sigma2 ≟ 0, S2fs0, S2f], If[sigma2 ≟ 0, S2nfs0, S2nf])
S1baseinv(x) = theta1 - (sigma1 * log(M1 / x - 1))
S2baseinv(x) = theta2 - (sigma2 * log(M2 / x - 1))
equilibria_nfnf: {Intersect[S1base, S2base, xmin, xmax]}
lim2top: S2base(theta2 + (6 * sigma2))
lim2bot: S2base(theta2 - (6 * sigma2))
lim1top: S1base(theta1 + (6 * sigma1))
lim1bot: S1base(theta1 - (6 * sigma1))
equilibria_nff = RemoveUndefined(Join[{{Intersect[S1base, S2baseinv, Min[lim2bot, lim2top], Max[lim2bot, lim2top]]}, If[(S1base(lim2top) > S2baseinv(lim2top)) ∧ (sigma2 > 0) ∨ (S1base(lim2top)  <  S2baseinv(lim2top)) ∧ (sigma2  <  0), {(M2, S1base(M2))}, {}], If[(S1base(lim2bot) > S2baseinv(lim2bot)) ∧ (sigma2  <  0) ∨ (S1base(lim2bot)  <  S2baseinv(lim2bot)) ∧ (sigma2 > 0), {(0, S1base(0))}, {}]}])
equilibria_fnf = RemoveUndefined(Join[{{Intersect[S2base, S1baseinv, Min[lim1bot, lim1top], Max[lim1bot, lim1top]]}, If[(S2base(lim1top) > S1baseinv(lim1top)) ∧ (sigma1 > 0) ∨ (S2base(lim1top)  <  S1baseinv(lim1top)) ∧ (sigma1  <  0), {(M1, S2base(M1))}, {}], If[(S2base(lim1bot) > S1baseinv(lim1bot)) ∧ (sigma1  <  0) ∨ (S2base(lim1bot)  <  S1baseinv(lim1bot)) ∧ (sigma1 > 0), {(0, S2base(0))}, {}]}])
equilibria_ff = Zip((y(a), x(a)), a, equilibria_nfnf)
equilibria = If(flip1 ≟ 1, If[sigma1 ≟ 0, If[flip2 ≟ 1, If[sigma2 ≟ 0, {}, {(S2base(theta1), theta1)}], If[sigma2 ≟ 0, {(theta2, theta1)}, {(S2baseinv(theta1), theta1)}]], If[flip2 ≟ 1, If[sigma2 ≟ 0, {(S1base(theta2), theta2)}, equilibria_ff], If[sigma2 ≟ 0, {(theta2, S1baseinv(theta2))}, equilibria_fnf]]], If[sigma1 ≟ 0, If[flip2 ≟ 1, If[sigma2 ≟ 0, {(theta1, theta2)}, {(theta1, S2baseinv(theta1))}], If[sigma2 ≟ 0, {}, {(theta1, S2base(theta1))}]], If[flip2 ≟ 1, If[sigma2 ≟ 0, {(S1baseinv(theta2), theta2)}, equilibria_nff], If[sigma2 ≟ 0, {(theta2, S1base(theta2))}, equilibria_nfnf]]])
E1_sub = (51.906179925080515, 94.5379971465656)
E1close = Max(Sequence[If[Length[E1_sub - Element[equilibria, i]]  <  sqrt(dx^(2) + dy^(2)) / 2, i, 0], i, 1, Length[equilibria]])
E1 = DynamicCoordinates(E1_sub, If[x(E1_sub) > xmax, xmax, If[x(E1_sub)  <  xmin, xmin, If[E1close > 0, x(Element[equilibria, E1close]), x(E1_sub)]]], If[y(E1_sub) > ymax, ymax, If[y(E1_sub)  <  ymin, ymin, If[E1close > 0, y(Element[equilibria, E1close]), y(E1_sub)]]])
E1stable = true
E1type = 1
E1label = If(E1type ≟ 1, "SN", If[E1type ≟ 2, "SS", If[E1type ≟ 3, "C", If[E1type ≟ 4, "US", If[E1type ≟ 5, "UN", "SD"]]]])
E2_sub = (1.9886051360058412, 3.315175685303399)
E2close = Max(Sequence[If[Length[E2_sub - Element[equilibria, i]]  <  sqrt(dx^(2) + dy^(2)) / 2, i, 0], i, 1, Length[equilibria]])
E2 = DynamicCoordinates(E2_sub, If[x(E2_sub) > xmax, xmax, If[x(E2_sub)  <  xmin, xmin, If[E2close > 0, x(Element[equilibria, E2close]), x(E2_sub)]]], If[y(E2_sub) > ymax, ymax, If[y(E2_sub)  <  ymin, ymin, If[E2close > 0, y(Element[equilibria, E2close]), y(E2_sub)]]])
E2type = 1
E2stable = true
E2label = If(E2type ≟ 1, "SN", If[E2type ≟ 2, "SS", If[E2type ≟ 3, "C", If[E2type ≟ 4, "US", If[E2type ≟ 5, "UN", "SD"]]]])
E3_sub = (35.21101749203677, 86.98347170937329)
E3close = Max(Sequence[If[Length[E3_sub - Element[equilibria, i]]  <  sqrt(dx^(2) + dy^(2)) / 2, i, 0], i, 1, Length[equilibria]])
E3 = DynamicCoordinates(E3_sub, If[x(E3_sub) > xmax, xmax, If[x(E3_sub)  <  xmin, xmin, If[E3close > 0, x(Element[equilibria, E3close]), x(E3_sub)]]], If[y(E3_sub) > ymax, ymax, If[y(E3_sub)  <  ymin, ymin, If[E3close > 0, y(Element[equilibria, E3close]), y(E3_sub)]]])
E3type = 5
E3stable = false
E3label = If(E3type ≟ 1, "SN", If[E3type ≟ 2, "SS", If[E3type ≟ 3, "C", If[E3type ≟ 4, "US", If[E3type ≟ 5, "UN", "SD"]]]])
n_e = 1 (slider from 0 to 3, step 1)
E1found: E1close > 0
E2found: (E1 ≠ E2) ∧ (E2close > 0)
E3found: (E1 ≠ E3) ∧ (E2 ≠ E3) ∧ (E3close > 0)
Etypes = First({E1type, E2type, E3type}, n_e)
xlabeldisplayed: "" + (LaTeX[xlabel]) + ""
ylabeldisplayed: "" + (LaTeX[ylabel]) + ""
null2label: "" + (LaTeX[ylabel]) + "\text{-nullcline}"
null1label: "" + (LaTeX[xlabel]) + "\text{-nullcline}"
gridx = Sequence(Segment[((i * dx), ymin), ((i * dx), ymax)], i, ceil(xmin / dx), floor(xmax / dx))
gridy = Sequence(Segment[(xmin, (i * dy)), (xmax, (i * dy))], i, ceil(ymin / dy), floor(ymax / dy))
dx2 = 50
dy2 = 50
gridx2 = Sequence(Segment[((i * dx2), ymin), ((i * dx2), ymax)], i, ceil(xmin / dx2), floor(xmax / dx2))
gridy2 = Sequence(Segment[(xmin, (i * dy2)), (xmax, (i * dy2))], i, ceil(ymin / dy2), floor(ymax / dy2))
Efound = Zip(If[a, 1, 0], a, First[{E1found, E2found, E3found}, n_e])
Estability = Zip(If[a, 1, 0], a, First[{E1stable, E2stable, E3stable}, n_e])