I am trying to plot a very basic bifurcation plot of my simple first order ODE.

using DifferentialEquations
using Gnuplot
using BifurcationKit
# this definition makes the dI a vector so have to send in vector initial conditions
# model = @ode_def begin
# dI = β*I*(1 - 1/𝕽 - I)
# end β 𝕽
odemodel(I, p, t) = p.β*I*(1 - 1/p.𝕽 - I)
prob = ODEProblem(odemodel, 1.0, (0.0, 500.0), (β = 0.04, 𝕽 = 1.2), jac=true)
sol = solve(prob)

I was told the example here would help me out but I am not sure how to get the jacobian of my system.

But I don’t have odefun… is that something I can create out of the ODEProblem object? For example, 𝕽 = 1.0 is a bifurcation point for a transcritical bifurcation and I was wondering if I could show that using BifurcationKit.jl