# ============================================================================
# Investigation of Allen's SI/SIR Epidemic Model Conjectures
# Focus: Equation (3) and Endemic Equilibrium Analysis
# 
# Paper: Ladas & Allen (2001) "Open Problems and Conjectures: SI and SIR 
#        Epidemic Models", J. Difference Equations and Applications, 7:5
# ============================================================================

# Step 1: Load DMB.txt library
read "DMB.txt";

# Step 2: Load plots package and set parameters (R0 > 1)
with(plots):
a1 := 1.0:  b1 := 1/3:  c1 := 1/3:
R0 := evalf(a1/(b1 + c1));

# Step 3: Define procedure to solve Equation (3) from the paper
# Equation (3a): (b+c)*x* = y*(1 - exp(-a*x*))
# Equation (3b): y* = 1 - x*(1 + c/b)
EndemicEq3 := proc(a, b, c)
    local eq1, xs, ys, sol:
    eq1 := (b + c)*xs = (1 - xs*(1 + c/b))*(1 - exp(-a*xs)):
    sol := fsolve(eq1, xs, 0.0001..0.9999):
    if sol = NULL then RETURN([0, 1]): fi:
    ys := 1 - sol*(1 + c/b):
    [evalf(sol), evalf(ys)]:
end:

# Step 4: Compute the endemic equilibrium from Equation (3)
endemic_pt := EndemicEq3(a1, b1, c1);

# Step 5: Verify Equation (3) is satisfied
xs := endemic_pt[1]:  ys := endemic_pt[2]:
LHS := (b1 + c1)*xs;
RHS := ys*(1 - exp(-a1*xs));

# Step 6: Define transformation and test SSg/SSSg from DMB.txt
T := AllenSIR(a1, b1, c1, x, y);
SSg(T, [x, y]);

# Step 7: Test SSSg (also fails due to transcendental equation)
SSSg(T, [x, y]);

# Step 8: Compute orbits from different initial conditions (Conjecture 1 test)
orb1 := ORB(T, [x, y], [0.1, 0.8], 0, 500):
orb2 := ORB(T, [x, y], [0.3, 0.5], 0, 500):
orb3 := ORB(T, [x, y], [0.05, 0.9], 0, 500):
orb1[-1]; orb2[-1]; orb3[-1];

# Step 9: Create phase portrait showing convergence to endemic equilibrium
p1 := pointplot(orb1, color = blue, symbol = point):
p2 := pointplot(orb2, color = red, symbol = point):
p3 := pointplot(orb3, color = green, symbol = point):
p_end := pointplot([endemic_pt], symbol = solidcircle, symbolsize = 15, color = black):
display(p1, p2, p3, p_end, title = "Convergence to Eq(3) Endemic Equilibrium");

# Step 10: Test SI model (c = 0 case - this is PROVEN in the paper)
T_SI := AllenSIR(0.8, 0.3, 0, x, y);
orb_SI := ORB(T_SI, [x, y], [0.1, 0.9], 0, 300):
orb_SI[-1];

# Step 11: Jacobian analysis for stability (using JAC from DMB.txt)
J := JAC(T, [x, y]);
J_endemic := subs({x = endemic_pt[1], y = endemic_pt[2]}, J);
Eigenvalues(Matrix(evalf(J_endemic)));

# Step 12: Compute eigenvalue magnitude (STABILITY PROOF)
# If |lambda| < 1, the endemic equilibrium is locally asymptotically stable
abs(0.7438152229 + 0.1956596179*I);

# ============================================================================
# RESULTS SUMMARY:
# ============================================================================
# Endemic equilibrium from Eq(3): (0.1423746737, 0.7152506526)
# Equation (3) verified: LHS = RHS = 0.09491644913
# All orbits converge to endemic equilibrium (Conjecture 1 evidence)
# Eigenvalues: 0.7438 +/- 0.1957i (complex conjugates)
# |lambda| = 0.7691 < 1 --> LOCALLY ASYMPTOTICALLY STABLE
# CONJECTURE 1 appears TRUE based on numerical investigation
# ============================================================================