Question about model change in Nonlinear system solve

I am using Julia to solve a system of nonlinear equations, but my model has phase changes, meaning there are switches in the equations of the model itself. If the solution enters a single-phase region, the original conservation equations should be changed to a form without the gas term, and if it is a two-phase region, it should be the normal conservation equations. What should I do?

My question is a steady-state solution problem.

One way to do that is a continuous relaxation of the discrete behavior

A small question, can this operation be set in a steady-state nonlinear problem? Or where can one find relevant examples?

function NonlinearSystem(eqs, unknowns, ps;
        observed = [],
        name = nothing,
        description = "",
        default_u0 = Dict(),
        default_p = Dict(),
        defaults = _merge(Dict(default_u0), Dict(default_p)),
        guesses = Dict(),
        initializesystem = nothing,
        initialization_eqs = Equation[],
        systems = NonlinearSystem[],
        connector_type = nothing,
        continuous_events = nothing, # this argument is only required for ODESystems, but is added here for the constructor to accept it without error
        discrete_events = nothing,   # this argument is only required for ODESystems, but is added here for the constructor to accept it without error
        checks = true,
        parameter_dependencies = Equation[],
        metadata = nothing,
        gui_metadata = nothing)
    continuous_events === nothing || isempty(continuous_events) ||
        throw(ArgumentError("NonlinearSystem does not accept `continuous_events`, you provided $continuous_events"))
    discrete_events === nothing || isempty(discrete_events) ||
        throw(ArgumentError("NonlinearSystem does not accept `discrete_events`, you provided $discrete_events"))
    name === nothing &&
        throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
    length(unique(nameof.(systems))) == length(systems) ||
        throw(ArgumentError("System names must be unique."))
    (isempty(default_u0) && isempty(default_p)) ||
        Base.depwarn(
            "`default_u0` and `default_p` are deprecated. Use `defaults` instead.",
            :NonlinearSystem, force = true)

    # Accept a single (scalar/vector) equation, but make array for consistent internal handling
    if !(eqs isa AbstractArray)
        eqs = [eqs]
    end

    # Copy equations to canonical form, but do not touch array expressions
    eqs = [wrap(eq.lhs) isa Symbolics.Arr ? eq : 0 ~ eq.rhs - eq.lhs for eq in eqs]

    jac = RefValue{Any}(EMPTY_JAC)

    ps′ = value.(ps)
    dvs′ = value.(unknowns)
    parameter_dependencies, ps′ = process_parameter_dependencies(
        parameter_dependencies, ps′)

    defaults = Dict{Any, Any}(todict(defaults))
    guesses = Dict{Any, Any}(todict(guesses))
    var_to_name = Dict()
    process_variables!(var_to_name, defaults, guesses, dvs′)
    process_variables!(var_to_name, defaults, guesses, ps′)
    process_variables!(
        var_to_name, defaults, guesses, [eq.lhs for eq in parameter_dependencies])
    process_variables!(
        var_to_name, defaults, guesses, [eq.rhs for eq in parameter_dependencies])
    defaults = Dict{Any, Any}(value(k) => value(v)
    for (k, v) in pairs(defaults) if v !== nothing)
    guesses = Dict{Any, Any}(value(k) => value(v)
    for (k, v) in pairs(guesses) if v !== nothing)

    isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))

    NonlinearSystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
        eqs, dvs′, ps′, var_to_name, observed, jac, name, description, systems, defaults,
        guesses, initializesystem, initialization_eqs, connector_type, parameter_dependencies,
        metadata, gui_metadata, checks = checks)
end

It’s not an option, you’d need to do it by hand for now, for example writing if as tanh