I have a set of nonlinear equations (which happen to be steady-state DEs) that I am trying to solve that have multiple real roots and 3 variables (x, y, & z in the code below)

In Mathematica I can get the multiple roots as I’d expect:

(the difference I care about in comparing the two graphs is that the curves go jagged though they shouldn’t. the axes are slightly different since I re-scaled in the Julia figure, and the points are data which can be neglected.)

I would like to find the multiple solutions I have in Mathematica. Naturally, JuMP finds the single optimal according to the specified objective, so this is as expected but not as desired.

I tried IntervalRootFinding.jl:

function core_eqns(u)
x, y, z = u
φ_R = 2 * 0.30 * x / (x + 0.8) + 0.05
φ_F = 0.35 - φ_R
φ_M = φ_F * 0.99
φ_AR = φ_F * 0.01
C_Rf = 63.7e-6 * φ_R - y
φ_Rf = φ_R - y / 63.7e-6
λ = 10 * x / (x + 0.03) * φ_Rf
# steady-state equations
A = 5 * φ_M - (1 + x) * λ
B = 2.25e6 * z * C_Rf - 3.78 * y - λ * y
C = 800 * (5e-4 - z) - 35 * 60^2 * φ_AR * 536e-6 * z / (z + 15e-6) + 3.78 * y - 2.25e6 * C_Rf * z - λ * z
return [A, B, C]
end
function SS_solver(p)
A_range = 0..15
B_range = 0..35e-6
C_range = 0..5e-4
interval = A_range × B_range × C_range
function f(x)
out = core_eqns(x, p)
return SVector(
out[1],
out[2],
out[3]
)
end
return roots(f, interval)
end

this gave me ~2000 intervals all with the flag :unknown, and it took 455s to run – which is quite undesirable since each the curves above are generated from computing the roots 100 times (and for comparison it takes <4s for all of the solutions to run in JuMP)

the DE steady solvers are the best combination of speed and robustness that I’ve found, but haven’t been able to get the multiple real solutions, I might just need to more systematically sweep different initial conditions

┌─ Curve type: EquilibriumCont
├─ Number of points: 2
├─ Type of vectors: Vector{Float64}
├─ Parameter p starts at 1.0e-5, ends at 0.001
├─ Algo: PALC
└─ Special points:
- # 1, endpoint at p ≈ +0.00100000,

HomotopyContinuation.jl is working for single solutions but I am still working out how to run the loop so I can plot the multiple solutions nicely like in the original question. Strictly speaking, HomotopyContinuation.jl has satisfied my question as originally stated, but I still need to get the analogous plot which requires looping through the solution multiple times (like is intended in the bifurcation analysis)