# HomotopyContinuation: Duffing Revisited: Plotting Solution from Multiple Paths

**URL:** <https://discourse.julialang.org/t/homotopycontinuation-duffing-revisited-plotting-solution-from-multiple-paths/124735>\
**Category:** Modelling & Simulations\
**Tags:** question, package\
**Created:** [January 13, 2025, 5:54pm UTC](https://discourse.julialang.org/t/homotopycontinuation-duffing-revisited-plotting-solution-from-multiple-paths/124735 "2025-01-13T17:54:18Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [January 13, 2025, 5:54pm UTC](https://discourse.julialang.org/t/homotopycontinuation-duffing-revisited-plotting-solution-from-multiple-paths/124735/1 "2025-01-13T17:54:18Z")

</div>

I am seeking to plot the solutions that HomotopyContinuation.jl returns in case of multiple solution branches.

Same question reformulated: I am seeking to reproduce [BifurcationKit: Duffing Equation: Time-Harmonic: Single Harmonic](https://discourse.julialang.org/t/bifurcationkit-duffing-equation-time-harmonic-single-harmonic/123937) with HomotopyContinuation.jl

Below are two MWEs. The first is linear and works fine. The second is non-linear and fails (likely to my erroneous post-processing). The correct solution is provided by HarmonicBalance.jl by redefining the plot function (details sadly missing).

First MWE

**1/ Define linear and non-linear system**

```julia
@var A B p

m=1; freq = .5; om0=2*π*freq; ga=0.1; stiffnlin=1e-3; F0 = 1.;

Flin = System(
    [
        m*(om0^2-p[1]^2)*A + ga*p[1]*B,
        -ga*p[1]*A + m*(om0^2-p[1]^2)*B - F0,        
    ];
    parameters = [p]
)

F = System(
    [
       m*(om0^2-p[1]^2)*A + ga*p[1]*B + stiffnlin*.75*(A^3 + A*B^2),
       -ga*p[1]*A + m*(om0^2-p[1]^2)*B + stiffnlin*.75*(B^3 + A^2*B) - F0,
        
    ];
    parameters = [p]
)

```

**2/ Solve and plot the linear system**

```julia
# solve system for multiple values of driving frequency using start_parameters and target_parameters
omdvec = Vector(0:0.01:10.)
Avec = zeros(length(omdvec))
Bvec = zeros(length(omdvec))
for (i,p) in enumerate(omdvec) 
    # Compute all solutions for F_p
    results = solve(Flin, target_parameters = [p])
    Avec[i] = real_solutions(results)[1][1]
    Bvec[i] = real_solutions(results)[1][2]
end

plot(omdvec,Avec)
plot!(omdvec,Bvec)
plot!(omdvec,sqrt.(Avec.^2+Bvec.^2))

```

Second MWE

**3/ Solve and plot the non-linear system**

Replace “Flin” by “F” in solve function above.

PS (for what it is worth): The page [How to solve a system of polynomial equations](https://www.juliahomotopycontinuation.org/guides/solve-first-system/) contains a tiny typo as realsolutions (without hyphen) should be real\_solutions (with hyphen) instead.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [January 14, 2025, 9:24pm UTC](https://discourse.julialang.org/t/homotopycontinuation-duffing-revisited-plotting-solution-from-multiple-paths/124735/2 "2025-01-14T21:24:11Z")

</div>

I did improved using

```julia
# solve system for multiple values of driving frequency using start_parameters and target_parameters
omdvec = Vector(0:0.01:10.)
Avec = Float64[]
Bvec = Float64[] 
start_solution = zeros(2)
start_p = 0. 
for (i,p) in enumerate(omdvec) 
    # Compute all solutions for F_p
    results = solve(F, start_solution; target_parameters = [p])
    if !(isempty(real_solutions(results))) 
        push!(Avec, real_solutions(results)[1][1]) 
        push!(Bvec, real_solutions(results)[1][2]) 
        start_solution = [real_solutions(results)[1][1], real_solutions(results)[1][2]]
        start_p = p  
    end 
end

plot(Avec)
plot!(Bvec)
plot!(sqrt.(Avec.^2+Bvec.^2))

```

to obtain

 ![Screenshot 2025-01-14 at 22.21.17](https://global.discourse-cdn.com/julialang/original/3X/8/b/8b870c70ca5eb79c16956c20bf04ac22453ca104.png)

I am, however, not quite there yet.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [January 14, 2025, 9:51pm UTC](https://discourse.julialang.org/t/homotopycontinuation-duffing-revisited-plotting-solution-from-multiple-paths/124735/3 "2025-01-14T21:51:49Z")

</div>

I miss functionality to plot multiple solutions for particular parameter values.

A solution is provided in HarmonicBalance.jl at [HarmonicBalance.jl/src/plotting\_Plots.jl at master · NonlinearOscillations/HarmonicBalance.jl · GitHub](https://github.com/NonlinearOscillations/HarmonicBalance.jl/blob/master/src/plotting_Plots.jl)

Are there other solutions? Thx.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [January 15, 2025, 5:03pm UTC](https://discourse.julialang.org/t/homotopycontinuation-duffing-revisited-plotting-solution-from-multiple-paths/124735/4 "2025-01-15T17:03:16Z")

</div>

Using

```julia
Avec = Vector{Float64}[]
Bvec = Vector{Float64}[]
Cvec = Vector{Float64}[]
for (i,ri) in enumerate(myrealsols)
    Aveci = []
    Bveci = [] 
    Cveci = [] 
    for (j,rj) in enumerate(ri)
        Avecj = rj[1]
        Bvecj = rj[2] 
        Cvecj = sqrt(Avecj^2+Bvecj^2)
        push!(Aveci,Avecj) 
        push!(Bveci,Bvecj) 
        push!(Cveci,Cvecj) 
    end 
    push!(Avec,Aveci)
    push!(Bvec,Bveci)
    push!(Cvec,Cveci)
end 

X = [i*ones(length(res)) for (i,res) in pairs(myrealsols)];

scatter(X, Cvec, ms=5, legend =:false)

```

I now obtain

 ![Screenshot 2025-01-15 at 17.58.21](https://global.discourse-cdn.com/julialang/original/3X/3/d/3d8a272ec160a251e43db14e7f1d496ae07398ac.png)

matching well with the previous solution.