second order linear response for branches without NaN

[oameye]

using HarmonicBalance
using Plots.Measures: mm
@variables α, ω, ω0, F, γ, t, x(t); # declare constant variables and a function x(t)

# define ODE
diff_eq = DifferentialEquation(d(x,t,2) + ω0*x + α*x^3 + γ*d(x,t) ~ F*cos(ω*t), x)

# specify the ansatz x = u(T) cos(ω*t) + v(T) sin(ω*t)
add_harmonic!(diff_eq, x, ω)

# implement ansatz to get harmonic equations
harmonic_eq = get_harmonic_equations(diff_eq)

fixed = (α => 1, ω0 => 1.0, γ => 0.005, F => 0.002)   # fixed parameters
varied = ω => range(0.95, 1.05, 100)           # range of parameter values
result = get_steady_states(harmonic_eq, varied, fixed)

plot_linear_response(result, x, 2, Ω_range=range(0.9,1.1,300), logscale=true, order=2)

BoundsError: attempt to access 300×37 Matrix{Float64} at index [1, 64]

Stacktrace:
 [1] setindex!
   @ .\array.jl:1024 [inlined]
 [2] get_linear_response(res::HarmonicBalance.Result{1, ComplexF64, Float64, FunctionWrappers.FunctionWrapper{Matrix{ComplexF64}, Tuple{Vector{ComplexF64}}}}, nat_var::Symbolics.Num, Ω_range::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, branch::Int64; order::Int64, show_progress::Bool)
   @ HarmonicBalance.LinearResponse C:\Users\User\.julia\packages\HarmonicBalance\GCsOJ\src\modules\LinearResponse\plotting.jl:112
 [3] get_linear_response
   @ C:\Users\User\.julia\packages\HarmonicBalance\GCsOJ\src\modules\LinearResponse\plotting.jl:89 [inlined]
 [4] plot_linear_response(res::HarmonicBalance.Result{1, ComplexF64, Float64, FunctionWrappers.FunctionWrapper{Matrix{ComplexF64}, Tuple{Vector{ComplexF64}}}}, nat_var::Symbolics.Num, branch::Int64; Ω_range::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, order::Int64, logscale::Bool, show_progress::Bool, kwargs::@Kwargs{})
   @ PlotsExt C:\Users\User\.julia\packages\HarmonicBalance\GCsOJ\ext\PlotsExt\linear_response.jl:37
 [5] top-level scope
   @ In[16]:1