# Brute force bifurcation diagram - need help with stiff solver

**URL:** <https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387>\
**Category:** Modelling & Simulations\
**Created:** [February 10, 2023, 12:42am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387 "2023-02-10T00:42:17Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [February 10, 2023, 12:42am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/1 "2023-02-10T00:42:17Z")

</div>

Hi All,

I’m trying to reproduce a bifurcation diagram (brute force) computed using XPP-AUT. My code is [here](https://github.com/epirecipes/chaos-measles/tree/main/krylova2013) and the original XPP code from Krylova and Earn (2013) is [here](https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2013.0098&file=rsif20130098supp2.pdf).

The XPP code seems to give a nice bifurcation plot:

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/8/0/80ee1f51cdfc34cf53bca5bde1e6650e0e16f922.png)

However, I’ve tried multiple solvers, and many of them fail with this very stiff problem. Here’s the best I’ve done so far, using a PositiveDomain callback:

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/2/b/2bbeff733169679c40c467d5bc3d17872597115e.png)

You can see gaps in the paths which aren’t present in the XPP code (which doesn’t throw any errors). Any idea of how to potentially improve the Julia code to be more robust?

To save you a trip to GitHub, here’s the Julia code:

```julia
using ModelingToolkit
using OrdinaryDiffEq
using DiffEqCallbacks
using DataFrames
using CSV

@parameters t R₀ γ μ a
@variables S(t) I(t) β(t)
D = Differential(t)

eqs = [D(S) ~ μ - β*S*I - μ*S,
       D(I) ~ β*S*I - (γ+μ)*I,
       β ~ R₀*(γ+μ)*(1+a*cos(2*π*t))]

@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)

p = [μ => 0.02, γ => 28.08, R₀ => 17.0, a => 0.08]
u₀ = [S => 0.9, I => 0.001]
tmin = 0.0
tmax = 650
transient = 600
strobe = 1.0
tspan = (tmin, tmax)
R0vec = collect(1.0:0.01:30.0)
solver = Tsit5()

prob = ODEProblem(simpsys, u₀, tspan, p; jac=true)

function prob_func(prob, i, repeat)
    remake(prob, p=[μ => 0.02, γ => 28.08, R₀ => R0vec[i], a => 0.08])
end

function output_func(sol, i)
    strobetimes = collect(transient:strobe:tmax)
    df = DataFrame(sol(strobetimes))
    rename!(df,[:t, :S, :I])
    df[!,"R0"] = fill(R0vec[i],length(strobetimes))
    df[!,"LogS"] = log10.(abs.(df.S))
    df[!,"LogI"] = log10.(abs.(df.I))
    return (df, false)
end

ensemble_prob = EnsembleProblem(prob,
                                prob_func = prob_func,
                                output_func = output_func)

@time sim = solve(ensemble_prob,
                  solver,
                  EnsembleThreads(),
                  trajectories = length(R0vec);
                  callback = PositiveDomain())

results = vcat(sim...)

CSV.write("bruteforceSIR_R0.dat", results; delim=" ", header=false)

using StatsPlots
using LaTeXStrings

p = @df results scatter(:R0,
                    :LogI,
                    xlabel=L"R_0",
                    ylabel=L"log_{10} I",
                    markersize=1.0,
                    color=:gray,
                    legend=false,
                    ylims=(-6,-2))
savefig(p, "bruteforceSIR_R0_statsplots.png")

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 10, 2023, 2:56am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/2 "2023-02-10T02:56:30Z")

</div>

if this is a stiff problem, why are you using tsit5? does Rodas5P do better?

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [February 10, 2023, 3:17am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/3 "2023-02-10T03:17:29Z")

</div>

I get lots of

```julia
Warning: Instability detected. Aborting

```

with Rodas5P (without changing anything else). Hence, I get a plot like this:

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/4/a/4ad77ba1c38f0a69adecfb3e10b4c523e5d61c44.png)

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 10, 2023, 3:22am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/4 "2023-02-10T03:22:47Z")

</div>

That makes very little sense. What about `Rosenbrock23` or `TRBDF2`?

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [February 10, 2023, 3:44am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/5 "2023-02-10T03:44:14Z")

</div>

They’re worse! Here’s TRBDF2 (again, using default settings as above). I agree it doesn’t make sense…

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/d/7/d733db18da65a23862484be9d56fd704d223d770.png)

Just for fun, here’s RK4

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/6/c/6c8838c0b3ec4628253de2cf47ae508731f4cb68.png)

Agreed that this doesn’t make a great deal of sense on first impression.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 10, 2023, 7:07am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/6 "2023-02-10T07:07:23Z")

</div>

Did you try lsoda or cvode\_bdf? Although i would be surprised if juliabdoes not beat them; i would expext xpp to use one of them

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 10, 2023, 7:59am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/7 "2023-02-10T07:59:35Z")

</div>

I would also say that you need multiple trials per `R0`

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [February 10, 2023, 8:30am UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/8 "2023-02-10T08:30:51Z")

</div>

I think the problem is the behavior of `PositiveDomain`. Its documentations says

> It reduces the next step by a certain scale factor until the extrapolated value at the next  
> time point is non-negative with a certain tolerance. Extrapolations are cheap to compute but might be inaccurate, so if a time step is changed it is additionally reduced by a safety factor of 0.9. Since extrapolated values are only non-negative up to a certain tolerance and in addition actual calculations might lead to negative values, also any negative values at the current time point are set to 0.

The last point is crucial. We want to see very small concentrations up to 1e-6 or so, but looking at the output data, they are simply set to zero by the callback.

Using `isoutofdomain=(u,p,t) -> any(x -> x < 0, u)` instead of `callback = PositiveDomain()` looks a lot better, filling out the “gaps” in the diagram. I also lowered the tolerances to `reltol=1e-11, abstol=1e-11` just to be sure, but that’s probably too low (and `isoutofdomain` should have a similar effect).

UPDATE: A combination of low `abstol` and the `isoutofdomain` check seems to be enough – no need for such a low `reltol` which just slows everything down.

![bruteforceSIR_R0_statsplots](https://global.discourse-cdn.com/julialang/original/3X/3/1/313a8c8fefe9dfbc4bb18fa0ee4e2e5de0f47e4e.png)

PS: Lowering the tolerance inside `PositiveDomain` might also have the same effect, didn’t try that yet.

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [February 10, 2023, 5:47pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/9 "2023-02-10T17:47:14Z")

</div>

> [@Sevi](#):
>
> eltol=1e-11, abstol=1e-11

Thanks! This looks a lot better (though not as good as the XPP code). I’ll dig in to the settings in XPP (I think it uses a Dormand-Prince scheme as default).

Any idea why ‘non-stiff’ solvers do better here? The system has an extreme transient at the beginning, which I can probably overcome with some fiddly initial condition wrangling (but this doesn’t seem necessary in XPP).

Here’s the code so far with a non-stiff solver (RK4)

```julia
using ModelingToolkit
using OrdinaryDiffEq
using DataFrames
using CSV

@parameters t R₀ γ μ a
@variables S(t) I(t) β(t)
D = Differential(t)

eqs = [D(S) ~ μ - β*S*I - μ*S,
       D(I) ~ β*S*I - (γ+μ)*I,
       β ~ R₀*(γ+μ)*(1+a*cos(2*π*t))]

@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)

p = [μ => 0.02, γ => 28.08, R₀ => 17.0, a => 0.08]
u₀ = [S => 0.9, I => 0.001]
tmin = 0.0
tmax = 650
transient = 600
strobe = 1.0
tspan = (tmin, tmax)
R0vec = collect(1.0:0.01:30.0)
solver = RK4()
tol = 1e-11

prob = ODEProblem(simpsys, u₀, tspan, p; jac=true) 

function prob_func(prob, i, repeat)
    remake(prob, p=[μ => 0.02, γ => 28.08, R₀ => R0vec[i], a => 0.08])
end

function output_func(sol, i)
    strobetimes = collect(transient:strobe:tmax)
    df = DataFrame(sol(strobetimes))
    rename!(df,[:t, :S, :I])
    df[!,"R0"] = fill(R0vec[i],length(strobetimes))
    df[!,"LogS"] = log10.(abs.(df.S))
    df[!,"LogI"] = log10.(abs.(df.I))
    return (df, false)
end

ensemble_prob = EnsembleProblem(prob,
                                prob_func = prob_func,
                                output_func = output_func)

@time sim = solve(ensemble_prob,
                  solver,
                  EnsembleThreads(),
                  trajectories = length(R0vec);
                  isoutofdomain=(u,p,t) -> any(x -> x < 0, u),
                  abstol=tol)

results = vcat(sim...)

CSV.write("bruteforceSIR_R0.dat", results; delim=" ", header=false)

using StatsPlots
using LaTeXStrings

p = @df results scatter(:R0,
                    :LogI,
                    xlabel=L"R_0",
                    ylabel=L"log_{10} I",
                    markersize=1.0,
                    color=:gray,
                    legend=false,
                    ylims=(-6,-2))
savefig(p, "bruteforceSIR_R0_statsplots.png")

```

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/2/e/2e95a0ac78c864fa9afd1fbc8e1cc1bb1c1d9689.png)

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [February 10, 2023, 8:07pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/10 "2023-02-10T20:07:13Z")

</div>

I did some more investigation, and I see what’s going on with the stiff solvers; if you set a=0, the system has two stable states, one with S=1 and I=0, and one with S\<1 and I\>0. Some of the stiff solvers seem to get stuck in the first, while others fail.

Here’s some example code with different solvers.

```julia
using ModelingToolkit
using OrdinaryDiffEq
using LSODA
using Sundials

@parameters t R₀ γ μ a
@variables S(t) I(t) β(t)
D = Differential(t)

eqs = [D(S) ~ μ - β*S*I - μ*S,
       D(I) ~ β*S*I - (γ+μ)*I,
       β ~ R₀*(γ+μ)*(1+a*cos(2*π*t))]

@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)

u₀ = [S => 0.9, I => 0.001]
tmin = 0.0
tmax = 650
tol = 1e-11
maxiters = 1e7

p₀ = [μ => 0.02, γ => 28.08, R₀ => 17.0, a => 0.0]

prob = ODEProblem(simpsys, u₀, (tmin, tmax), p₀; jac=true)
sol_rk4 = solve(prob, RK4();
                maxiters = maxiters,
                isoutofdomain=(u,p,t) -> any(x -> (x <= 0 || x >= 1), u),
                abstol=tol)

sol_lsoda = solve(prob, lsoda();
                maxiters = maxiters,
                abstol=tol)

sol_rosenbrock23 = solve(prob, Rosenbrock23();
                maxiters = maxiters,
                isoutofdomain=(u,p,t) -> any(x -> (x <= 0 || x >= 1), u),
                abstol = tol)

sol_cvodebdf = solve(prob, CVODE_BDF())

```

Here’s the output of the final state.

```julia
julia> print(sol_rk4.u[end])
[0.05882344413440804, 0.0006699054322332516]
julia> print(sol_lsoda.u[end])
[0.05882352941176473, 0.0006698764915218742]
julia> print(sol_rosenbrock23.u[end])
[0.9999994236820176, 1.604027727172643e-47]

```

CVODE\_BDF is interesting. It works if I don’t set the tolerance to a low value.

```julia
julia> sol_cvodebdf = solve(prob, CVODE_BDF())
retcode: Success
Interpolation: 3rd order Hermite
t: 2679-element Vector{Float64}:
   0.0
   9.365055056001597e-5
   0.00018730110112003193
   0.00045236993728053627
   ⋮
 649.6776241741102
 649.912024827687
 650.0
u: 2679-element Vector{Vector{Float64}}:
 [0.9, 0.001]
 [0.8999583517992265, 0.001039101101715838]
 [0.8999150700341084, 0.0010797289334348362]
 [0.8997831768370831, 0.0012034963123768742]
 ⋮
 [0.05857054135397698, 0.0006870603487333439]
 [0.05853834794611831, 0.0006657353926255108]
 [0.05856604433672487, 0.0006581134035398938]

```

but fails if I drop the tolerance.

```julia
julia> sol_cvodebdf = solve(prob, CVODE_BDF(); abstol=1e-11)

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 2.9051e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 2.9051e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 2.9051e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 2.9051e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 1.79233e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 1.79233e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 1.79233e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 1.79233e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 1.11224e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  Internal t = 5.43326 and h = 1.11224e-16 are such that t + h = t on the next step. The solver will continue anyway.

[CVODES WARNING] CVode
  The above warning has been issued mxhnil times and will not be issued again for this problem.

retcode: MaxIters
Interpolation: 3rd order Hermite
t: 100001-element Vector{Float64}:
 0.0
 6.622041520638907e-5
 0.00013244083041277815
 0.0003449609238189219
 ⋮
 5.433256962605888
 5.433256962605888
 5.433256962605888
u: 100001-element Vector{Vector{Float64}}:
 [0.9, 0.001]
 [0.8999708850285875, 0.001027335790869192]
 [0.8999409716210338, 0.001055417802667344]
 [0.899839159128681, 0.0011509726910748413]
 ⋮
 [6.134513963058041e152, -6.13451396305804e152]
 [6.134513963058041e152, -6.13451396305804e152]

```

If I change R0 to 2.0 (which should make the system _less_ stiff), then lsoda and CVODE\_BDF fails too:

```julia
julia> print(sol_rk4.u[end])
[0.5000006296727381, 0.00035587073771302234]
julia> print(sol_lsoda.u[end])
[NaN, NaN]
julia> print(sol_rosenbrock23.u[end])
[0.9999995646381588, 5.974907074363282e-99]
julia> print(sol_cvodebdf.u[end])
[1.7885027908424116e153, -1.788502790842412e153]

```

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [February 10, 2023, 10:41pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/11 "2023-02-10T22:41:23Z")

</div>

I still don’t know why the stiff solvers can’t handle the transient properly, but with initial conditions closer to the equilibrium cycle, stiff solvers like `Rosenbrock23` now give better results (though why it’s not working still for low R0 is odd to me). Code below.

```julia
using ModelingToolkit
using OrdinaryDiffEq
using DataFrames
using CSV

@parameters t R₀ γ μ a
@variables S(t) I(t) β(t)
D = Differential(t)

eqs = [D(S) ~ μ - β*S*I - μ*S,
       D(I) ~ β*S*I - (γ+μ)*I,
       β ~ R₀*(γ+μ)*(1+a*cos(2*π*t))]

@named sys = ODESystem(eqs)
simpsys = structural_simplify(sys)

S₀ = 1/R₀
I₀ = (μ/(μ+γ))*(1-S₀)
u₀ = [S => S₀, I => I₀]
p = [μ => 0.02, γ => 28.08, R₀ => 17.0, a => 0.08]

tmin = 0.0
tmax = 650
transient = 600
strobe = 1.0
tspan = (tmin, tmax)
R0vec = collect(1.0:0.01:30.0)
solver = Rosenbrock23()
tol = 1e-11
maxiters = 1e7

prob = ODEProblem(simpsys, u₀, tspan, p; jac=true) 

# Utility functions

function prob_func(prob, i, repeat)
    remake(prob, p=[μ => 0.02, γ => 28.08, R₀ => R0vec[i], a => 0.08])
end

function output_func(sol, i)
    strobetimes = collect(transient:strobe:tmax)
    df = DataFrame(sol(strobetimes))
    rename!(df,[:t, :S, :I])
    df[!,"R0"] = fill(R0vec[i],length(strobetimes))
    df[!,"LogS"] = log10.(abs.(df.S))
    df[!,"LogI"] = log10.(abs.(df.I))
    return (df, false)
end

# Run ensemble

ensemble_prob = EnsembleProblem(prob,
                                prob_func = prob_func,
                                output_func = output_func)

@time sim = solve(ensemble_prob,
                  solver,
                  EnsembleThreads(),
                  trajectories = length(R0vec);
                  maxiters = maxiters,
                  isoutofdomain=(u,p,t) -> any(x -> x < 0, u),
                  abstol=tol)

# Output

results = vcat(sim...)
CSV.write("bruteforceSIR_R0.dat", results; delim=" ", header=false)

# Plotting

using StatsPlots
using LaTeXStrings

p = @df results scatter(:R0,
                    :LogI,
                    xlabel=L"R_0",
                    ylabel=L"log_{10} I",
                    markersize=1.0,
                    color=:gray,
                    legend=false,
                    ylims=(-6,-2))
savefig(p, "bruteforceSIR_R0_statsplots.png")

```

 ![bruteforceSIR_R0](https://global.discourse-cdn.com/julialang/original/3X/1/9/19f6c0e4100d06834e8e22cae4838845108badf9.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 11, 2023, 1:43pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/12 "2023-02-11T13:43:59Z")

</div>

> [@sdwfrost](#):
>
> I still don’t know why the stiff solvers can’t handle the transient properly

They are not stability bound but tolerance bound, so they might be taking large steps where the difference between two trajectories is within the tolerances. RK methods are forced to take really small steps and thus solve below the required error tolerance.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [August 22, 2023, 2:21pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/13 "2023-08-22T14:21:02Z")

</div>

I found a trick for a different context. Use I=exp(Ilog) and

```julia
function sir!(du, X, p, t = 0)
    @unpack R₀, γ, μ, a = p
    S, Ilog, u, v = X
    I = exp(Ilog)
    β = R₀ * (γ + μ) * (1 + a * cos(2pi*t)
    du[1] = μ - β * S * I - μ * S
    du[2] = β * S - (γ + μ)
    du
end

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [July 12, 2024, 4:56pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/14 "2024-07-12T16:56:10Z")

</div>

@Datseris could be an interesting problem for your continuation procedure, no?

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [July 12, 2024, 9:47pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/15 "2024-07-12T21:47:48Z")

</div>

yeah, from the original picture it appears that you only want the stable branches anyways. Just plug your system through our `global_continuation`. Use `AttractorsViaRecurrences` for the mapper, it doesn’t care if your problem is stiff or not; if the ODEsolve options you have chosen can evolve it, that is all that matters. See here if you are unfamiliar with `global_continuation`: [Attractors.jl Tutorial · Attractors.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/tutorial/)

note, with our method you will get the bifurcation diagram for both attractors that you mentioned exist in the system.

p.s.: It isn’t “my” procedure, it is “ours” 😉

* * *

Actually, you don’t need `global_continuation`. You just need to apply the `orbitdiagram` function. How do you collect points to plot the ODE solution in the first picture you showed? Is it a poincare crossing, a stroboscopic map? Or you just literally sfcatter plot the last X points of the timeseries? Here is `orbitdiagram`: [Orbit diagrams · ChaosTools.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/orbitdiagram/#ChaosTools.orbitdiagram) see the second example in this linked page that does it for a continuous time system. I think orbit diagram is what you call “brute force bifurcation diagram”.

---

<div class="post-metadata">

**Author:** ![magister-ludi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/magister-ludi/32/4003_2.png) [@magister-ludi](https://discourse.julialang.org/u/magister-ludi)\
**Post date:** [July 13, 2024, 1:33pm UTC](https://discourse.julialang.org/t/brute-force-bifurcation-diagram-need-help-with-stiff-solver/94387/16 "2024-07-13T13:33:46Z")

</div>

> [@Datseris](#):
>
> It isn’t “my” procedure, it is “ours”

👍