Merge 6a752f4 into 88675ff

docs/mkdocs.yml

@@ -27,3 +27,4 @@ pages:
     - examples.md
     - examples/calcium.md
     - examples/van_der_pol.md
+    - examples/bazykin_85.md

docs/src/examples.md

@@ -4,5 +4,6 @@
 Pages = [
     "examples/calcium.md",
     "examples/van_der_pol.md",
+    "examples/bazykin_85.md",
 ]
 ```

docs/src/examples/bazykin_85.md

@@ -0,0 +1,119 @@
+# Bazykin's predator-prey system
+
+```@example bazykin85
+using Bifurcations
+using Bifurcations.Examples: Bazykin85
+
+solver = init(Bazykin85.prob)
+solve!(solver)
+```
+
+Plot codimension-one bifurcations:
+
+```@example bazykin85
+using Bifurcations: plot, plot!  # workaround
+using Plots
+
+plt1 = plot(solver.sol)
+savefig(plt1, "bazykin85-1.png"); nothing # hide
+```
+
+![](bazykin85-1.png)
+
+Let's follow the Hopf and Saddle-Node bifurcations:
+
+```@example bazykin85
+using Bifurcations: Codim1, special_points
+using Setfield: @lens
+
+sn_point, = sort!(
+    special_points(solver, Codim1.PointTypes.saddle_node),
+    by=p->p.u0[1])
+hopf_point, = special_points(solver, Codim1.PointTypes.hopf)
+
+point_list = [sn_point, hopf_point]
+
+codim2_solvers = []
+for point in point_list
+    @show point
+
+    codim2_prob = BifurcationProblem(
+        point,
+        solver,
+        (@lens _.δ),
+        (0.0, 10.0),
+    )
+    codim2_solver = init(
+        codim2_prob;
+        nominal_angle_rad = 0.01,
+        max_samples = 1000,
+    )
+    push!(codim2_solvers, codim2_solver)
+    solve!(codim2_solver)
+
+    @show codim2_solver
+end
+
+sn_solver1, hopf_solver1 = codim2_solvers
+
+nothing # hide
+```
+
+Switch to Hopf bifurcation via Bogdanov-Takens bifurcation:
+
+```@example bazykin85
+using Bifurcations: Codim2
+bt_point, = sort(
+    special_points(sn_solver1,
+                   Codim2.PointTypes.bogdanov_takens);
+    by = p -> p.u0[end - 1],  # α
+    rev = true)
+hopf_prob = BifurcationProblem(bt_point, sn_solver1)
+hopf_solver2 = init(hopf_prob)
+solve!(hopf_solver2)
+@show hopf_solver2
+
+push!(codim2_solvers, hopf_solver2)
+
+nothing # hide
+```
+
+Switch to the fold bifurcation of limit cycle via Bautin bifurcation:
+
+```@example bazykin85
+using Bifurcations: Codim2
+using Bifurcations.Codim2LimitCycle: FoldLimitCycleProblem
+bautin_point_list = special_points(hopf_solver2, Codim2.PointTypes.bautin)
+@assert length(bautin_point_list) == 1
+bautin_point, = bautin_point_list
+flc_prob = FoldLimitCycleProblem(
+    bautin_point,
+    hopf_solver2;
+    num_mesh = 20,
+    degree = 3,
+)
+flc_solver = init(
+    flc_prob;
+    start_from_nearest_root = true,
+    max_branches = 0,
+)
+solve!(flc_solver)
+@show flc_solver
+
+push!(codim2_solvers, flc_solver)
+
+nothing # hide
+```
+
+Merge continuations and draw the bifurcation diagram:
+
+```@example bazykin85
+plt2 = plot()
+for s in codim2_solvers
+    n = length(s.sol.sweeps[1].super.u[1])
+    plot!(plt2, s, vars=(n - 1, n))
+end
+savefig(plt2, "bazykin85-2.png"); nothing # hide
+```
+
+![](bazykin85-2.png)

src/plotting.jl

@@ -438,7 +438,8 @@ function plot!(plt,
                                 BifurcationSolution,
                                 BifurcationSolver},
                args...;
-               resolve_points = plottable isa Codim1Solver,
+               resolve_points = plottable isa Union{Codim1Solver,
+                                                    Codim2Solver},
                include_points = !(plottable isa Codim1LCSolver),  # TODO: don't
                vars = length(args) > 0 ? nothing : default_vars(plottable),
                bif_style = STYLE,