How can I solve this 4-equation NonlinearSystem? It gives `retcode: Unstable`

a revised version below:

using NonlinearSolve
using ModelingToolkit

# !!! used UNITS !!!
#  eq/liter
#  Pa

# constants
mmHG = 133.322 # Pa
Kc = 2.46e-11/mmHG
Ka = 3e-7
@variables HA A H HCO3
@parameters pCO2 Atot SID

eqs = [
	H * A ~ Ka * HA
	H * HCO3 ~ Kc * pCO2
	Atot ~ HA + A
	SID ~ HCO3 + A
]

@mtkbuild ns_stewart = NonlinearSystem(eqs)
display(ns_stewart)
println()

# Initial values (guesses)
guesses_stewart = (
	HCO3 => 24.9e-3,
	H => 3.95e-8,
	A => 17.1e-3,
	HA => 2e-3
)

ps = (
    Atot => 19e-3,
    SID => 41.7e-3,
    pCO2 => 5.3e3
)

prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)
sol = solve(prob)

@show sol
@show H = -log10(sol[1])
println("HCO3 = $(sol[2]*1000) meq/liter")

NEW Results:

julia> include("stewart1.jl")
Model ns_stewart with 2 equations
Unknowns (2):
  H
  HCO3
Parameters (3):
  pCO2
  Atot
  SID

sol = [3.9275240698625016e-8, 0.02489948137605699]
H = -(log10(sol[1])) = 7.405881144702646
HCO3 = 24.89948137605699 meq/liter
1 Like