diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl
index 4a1e8fb1..14dc78fb 100644
--- a/src/HarmonicBalance.jl
+++ b/src/HarmonicBalance.jl
@@ -17,7 +17,6 @@ module HarmonicBalance
     export @variables
     export d
 
-
     import Base: ComplexF64, Float64; export ComplexF64, Float64
     ComplexF64(x::Complex{Num}) = ComplexF64(Float64(x.re) + im*Float64(x.im))
     Float64(x::Complex{Num}) = Float64(ComplexF64(x))
@@ -52,17 +51,18 @@ module HarmonicBalance
     include("saving.jl")
     include("transform_solutions.jl")
     include("plotting_Plots.jl")
-    include("hysteresis_sweep.jl")
 
     include("modules/HC_wrapper.jl")
     using .HC_wrapper
 
     include("modules/LinearResponse.jl")
     using .LinearResponse
+    export plot_linear_response, plot_rotframe_jacobian_response
 
     include("modules/TimeEvolution.jl")
     using .TimeEvolution
     export ParameterSweep, ODEProblem, solve
+    export plot_1D_solutions_branch, follow_branch
 
     include("modules/LimitCycles.jl")
     using .LimitCycles
diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl
index 9d105f41..a8094237 100644
--- a/src/HarmonicEquation.jl
+++ b/src/HarmonicEquation.jl
@@ -106,12 +106,16 @@ $(TYPEDSIGNATURES)
 Check if `eom` is rearranged to the standard form, such that the derivatives of the variables are on one side.
 """
 function is_rearranged(eom::HarmonicEquation)
-    tvar = get_independent_variables(eom)[1]
+    tvar = get_independent_variables(eom)[1]; dvar = d(get_variables(eom), tvar);
+    lhs = getfield.(eom.equations, :lhs); rhs = getfield.(eom.equations, :rhs);
 
-    # Hopf-containing equations are arranged by construstion (impossible to check)
-    isequal(getfield.(eom.equations, :rhs), d(get_variables(eom), tvar)) || in("Hopf", getfield.(eom.variables, :type))
-end
+    HB_bool = isequal(rhs, dvar)
+    hopf_bool = in("Hopf", getfield.(eom.variables, :type))
+    MF_bool = !any([occursin(str1,str2) for str1 in string.(dvar) for str2 in string.(lhs)])
 
+    # Hopf-containing equations or MF equation are arranged by construstion
+    HB_bool || hopf_bool || MF_bool
+end
 
 """
 $(TYPEDSIGNATURES)
diff --git a/src/hysteresis_sweep.jl b/src/hysteresis_sweep.jl
deleted file mode 100644
index 3078a308..00000000
--- a/src/hysteresis_sweep.jl
+++ /dev/null
@@ -1,101 +0,0 @@
-export plot_1D_solutions_branch
-
-
-"""Calculate distance between a given state and a stable branch"""
-function _closest_branch_index(res::Result,state::Vector{Float64},index::Int64)
-    #search only among stable solutions
-    physical     = filter_solutions.(
-                            filter_solutions.(
-                            res.solutions, res.classes["physical"]),res.classes["stable"])
-
-    steadystates = reduce(hcat,physical[index]) #change structure for easier indexing
-    distances    = vec(sum(abs2.(steadystates .- state),dims=1))
-    return argmin(replace!(distances, NaN=>Inf))
-end
-
-
-"""Return the indexes and values following stable branches along a 1D sweep.
-When a no stable solutions are found (e.g. in a bifurcation), the next stable solution is calculated by time evolving the previous solution (quench).
-
-Keyword arguments
-- `y`:  Dependent variable expression (parsed into Symbolics.jl) to evaluate the followed solution branches on .
-- `sweep`: Direction for the sweeping of solutions. A `right` (`left`) sweep proceeds from the first (last) solution, ordered as the sweeping parameter.
-- `tf`: time to reach steady
-- `ϵ`: small random perturbation applied to quenched solution, in a bifurcation in order to favour convergence in cases where multiple solutions are identically accessible (e.g. symmetry breaking into two equal amplitude states)
-"""
-function follow_branch(starting_branch::Int64,res::Result; y="sqrt(u1^2+v1^2)",sweep="right",tf=10000,ϵ=1E-4)
-    sweep_directions = ["left", "right"]
-    sweep ∈ sweep_directions  || error("Only the following (1D) sweeping directions are allowed:  ", sweep_directions)
-
-    #filter stable solutions
-    Y  = transform_solutions(res, y)
-    Ys = filter_solutions.(filter_solutions.(Y, res.classes["physical"]),res.classes["stable"]);
-    Ys = real.(reduce(hcat,Ys)) #facilitate indexing. Ys is an array (length(sweep))x(#solutions)
-
-    followed_branch    = zeros(Int64,length(Y)) #followed branch indexes
-    followed_branch[1] = starting_branch
-
-    if sweep == "left"
-        Ys = reverse(Ys,dims=2)
-    end
-
-    p1 = collect(keys(res.swept_parameters))[1] #parameter values
-
-    for i in 2:size(Ys)[2]
-        s = Ys[followed_branch[i-1],i] #solution amplitude in the current branch and current parameter index
-        if !isnan(s) #the solution is not unstable or unphysical
-            followed_branch[i] = followed_branch[i-1]
-        else #bifurcation found
-            if sweep == "right"
-                next_index = i
-            elseif sweep == "left" #reverse order of iterators
-                next_index = length(Y)-i+1
-            end
-
-            # create a synthetic starting point out of an unphysical solution: quench and time evolve
-            #the actual solution is complex there, i.e. non physical. Take real part for the quench.
-            sol_dict = get_single_solution(res, branch=followed_branch[i-1], index= next_index)
-            print("bifurcation @ ", p1 ," = ", real(sol_dict[p1])," ")
-            sol_dict   = Dict(zip(keys(sol_dict),real.(values(sol_dict)) .+ 0.0im .+ ϵ*rand(length(values(sol_dict)))))
-            problem_t  = TimeEvolution.ODEProblem(res.problem.eom, steady_solution=sol_dict, timespan=(0,tf))
-            res_t      = TimeEvolution.solve(problem_t,saveat=tf)
-
-            followed_branch[i] = _closest_branch_index(res,res_t.u[end],next_index) #closest branch to final state
-
-            print("switched branch ", followed_branch[i-1] ," -> ", followed_branch[i],"\n")
-        end
-    end
-
-    if sweep == "left"
-        Ys = reverse(Ys,dims=2)
-        followed_branch = reverse(followed_branch)
-    end
-
-    return followed_branch,Ys
-end
-
-
-"""1D plot with the followed branch highlighted"""
-function plot_1D_solutions_branch(starting_branch::Int64,res::Result; x::String, y::String, sweep="right",tf=10000,ϵ=1E-4, xscale="linear",yscale="linear",plot_only=["physical"],marker_classification="stable",filename=nothing,kwargs...)
-    plot_1D_solutions(res; x=x, y=y, xscale=xscale,yscale=yscale,
-                        plot_only=plot_only,marker_classification=marker_classification,filename=filename,kwargs...)
-
-    followed_branch,Ys = follow_branch(starting_branch,res,y=y,sweep=sweep,tf=tf,ϵ=ϵ)
-    if sweep == "left"
-        m = "<"
-    elseif sweep == "right"
-        m = ">"
-    end
-
-    Y_followed = [Ys[branch,param_idx] for (param_idx,branch) in enumerate(followed_branch)]
-
-    lines = plot(transform_solutions(res, x),Y_followed,c="k",marker=m; _set_Plots_default..., kwargs...)
-
-    #extract plotted data and return it
-    xdata,ydata = [line.get_xdata() for line in lines], [line.get_ydata() for line in lines]
-    markers = [line.get_marker() for line in lines]
-    save_dict= Dict(string(x) => xdata,string(y)=>ydata,"markers"=>markers)
-
-    !isnothing(filename) ? JLD2.save(_jld2_name(filename), save_dict) : nothing
-    return save_dict
-end
\ No newline at end of file
diff --git a/src/modules/LinearResponse/Lorentzian_spectrum.jl b/src/modules/LinearResponse/Lorentzian_spectrum.jl
index 8d66cb27..5bcb91b1 100644
--- a/src/modules/LinearResponse/Lorentzian_spectrum.jl
+++ b/src/modules/LinearResponse/Lorentzian_spectrum.jl
@@ -67,21 +67,23 @@ _get_as(hvars::Vector{HarmonicVariable}) = findall(x -> isequal(x.type, "a"), hv
 
 
 #   Returns the spectra of all variables in `res` for `index` of `branch`.
-function JacobianSpectrum(res::Result; index::Int, branch::Int)
+function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false)
+    hvars = res.problem.eom.variables # fetch the vector of HarmonicVariable
+    # blank JacobianSpectrum for each variable
+    all_spectra = Dict{Num, JacobianSpectrum}([[nvar, JacobianSpectrum([])] for nvar in getfield.(hvars, :natural_variable)])
 
-    res.classes["stable"][index][branch] || error("\nThe solution is unstable - it has no JacobianSpectrum!\n")
+    if force
+        res.classes["stable"][index][branch] || return all_spectra # if the solution is unstable, return empty spectra
+    else
+        res.classes["stable"][index][branch] || error("\nThe solution is unstable - it has no JacobianSpectrum!\n")
+    end
 
     solution_dict = get_single_solution(res, branch=branch, index=index)
-
-    hvars = res.problem.eom.variables # fetch the vector of HarmonicVariable
     λs, vs = eigen(res.jacobian(solution_dict))
 
     # convert OrderedDict to Dict - see Symbolics issue #601
     solution_dict = Dict(get_single_solution(res, index=index, branch=branch))
 
-    # blank JacobianSpectrum for each variable
-    all_spectra = Dict{Num, JacobianSpectrum}([[nvar, JacobianSpectrum([])] for nvar in getfield.(hvars, :natural_variable)])
-
     for (j, λ) in enumerate(λs)
         eigvec = vs[:, j] # the eigenvector
 
@@ -130,4 +132,3 @@ function _simplify_spectra!(spectra::Dict{Num, JacobianSpectrum})
         spectra[var] = _simplify_spectrum(spectra[var])
     end
 end
-
diff --git a/src/modules/LinearResponse/plotting.jl b/src/modules/LinearResponse/plotting.jl
index 877a48c7..2408d59b 100644
--- a/src/modules/LinearResponse/plotting.jl
+++ b/src/modules/LinearResponse/plotting.jl
@@ -1,9 +1,6 @@
 using Plots, Latexify, ProgressMeter
-using Latexify.LaTeXStrings
-import ..HarmonicBalance: dim, _get_mask
-using Plots.PlotUtils.ColorSchemes
-
-export plot_linear_response
+using HarmonicBalance: _set_Plots_default
+export plot_linear_response, plot_rotframe_jacobian_response
 
 
 function get_jacobian_response(res::Result, nat_var::Num, Ω_range, branch::Int; show_progress=true)
@@ -23,6 +20,21 @@ function get_jacobian_response(res::Result, nat_var::Num, Ω_range, branch::Int;
     end
     C
 end
+function get_jacobian_response(res::Result, nat_var::Num, Ω_range, followed_branches::Vector{Int}; show_progress=true, force=false)
+
+    spectra = [JacobianSpectrum(res, branch=branch, index = i, force=force) for (i, branch) in pairs(followed_branches)]
+    C = Array{Float64, 2}(undef,  length(Ω_range), length(spectra))
+
+    if show_progress
+        bar = Progress(length(CartesianIndices(C)), 1, "Diagonalizing the Jacobian for each solution ... ", 50)
+    end
+    # evaluate the Jacobians for the different values of noise frequency Ω
+    for ij in CartesianIndices(C)
+        C[ij] = abs(evaluate(spectra[ij[2]][nat_var], Ω_range[ij[1]]))
+        show_progress ? next!(bar) : nothing
+    end
+    C
+end
 
 
 function get_linear_response(res::Result, nat_var::Num, Ω_range, branch::Int; order, show_progress=true)
@@ -95,8 +107,23 @@ X = Vector{Float64}(collect(values(res.swept_parameters))[1][stable])
 C = order == 1 ? get_jacobian_response(res, nat_var, Ω_range, branch, show_progress=show_progress) : get_linear_response(res, nat_var, Ω_range, branch; order=order, show_progress=show_progress)
 C = logscale ? log.(C) : C
 
-heatmap(X, Ω_range,  C; color=:viridis,
-    xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
+xlabel = latexify(string(first(keys(res.swept_parameters)))); ylabel = latexify("Ω");
+heatmap(X, Ω_range,  C; color=:viridis, xlabel=xlabel, ylabel=ylabel, _set_Plots_default..., kwargs...)
+end
+function plot_linear_response(res::Result, nat_var::Num, followed_branches::Vector{Int}; Ω_range, logscale=false, show_progress=true, switch_axis = false, force=true, kwargs...)
+    length(size(res.solutions)) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
+
+    X = Vector{Float64}(collect(first(values(res.swept_parameters))))
+
+    C = get_jacobian_response(res, nat_var, Ω_range, followed_branches; force=force)
+    C = logscale ? log.(C) : C
+
+    xlabel = latexify(string(first(keys(res.swept_parameters)))); ylabel = latexify("Ω");
+    if switch_axis
+        heatmap(Ω_range, X, C'; color=:viridis, xlabel=ylabel, ylabel=xlabel, _set_Plots_default..., kwargs...)
+    else
+        heatmap(X, Ω_range, C; color=:viridis, xlabel=xlabel, ylabel=ylabel, _set_Plots_default..., kwargs...)
+    end
 end
 
 """
@@ -123,7 +150,7 @@ function plot_rotframe_jacobian_response(res::Result; Ω_range, branch::Int, log
     C = logscale ? log.(C) : C
 
     heatmap(X, Ω_range,  C; color=:viridis,
-       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
+       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), _set_Plots_default..., kwargs...)
 end
 
 function plot_eigenvalues(res; branch, type=:imag, projection= v -> 1, cscheme = :default, kwargs...)
diff --git a/src/modules/LinearResponse/response.jl b/src/modules/LinearResponse/response.jl
index 5552cd9c..c9a52d9c 100644
--- a/src/modules/LinearResponse/response.jl
+++ b/src/modules/LinearResponse/response.jl
@@ -1,7 +1,4 @@
 export get_response
-export plot_response
-
-
 
 """
     get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)
@@ -95,8 +92,3 @@ end
 # rotating frame into the lab frame
 _plusamp(uv) = norm(uv)^2 - 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
 _minusamp(uv) = norm(uv)^2 + 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
-
-
-
-
-
diff --git a/src/modules/TimeEvolution.jl b/src/modules/TimeEvolution.jl
index d543748b..6264601d 100644
--- a/src/modules/TimeEvolution.jl
+++ b/src/modules/TimeEvolution.jl
@@ -2,7 +2,9 @@ module TimeEvolution
 
     using ..HarmonicBalance
     using Symbolics
+    using Plots
     using OrdinaryDiffEq
+    using OrderedCollections
     using DSP
     using FFTW
     using Peaks
@@ -12,5 +14,6 @@ module TimeEvolution
     include("TimeEvolution/ODEProblem.jl")
     include("TimeEvolution/FFT_analysis.jl")
     include("TimeEvolution/sweeps.jl")
+    include("TimeEvolution/hysteresis_sweep.jl")
 
 end
diff --git a/src/modules/TimeEvolution/hysteresis_sweep.jl b/src/modules/TimeEvolution/hysteresis_sweep.jl
new file mode 100644
index 00000000..cb25035f
--- /dev/null
+++ b/src/modules/TimeEvolution/hysteresis_sweep.jl
@@ -0,0 +1,80 @@
+export plot_1D_solutions_branch, follow_branch
+
+"""Calculate distance between a given state and a stable branch"""
+function _closest_branch_index(res::Result, state::Vector{Float64}, index::Int64)
+    #search only among stable solutions
+    stable = HarmonicBalance._apply_mask(res.solutions, HarmonicBalance._get_mask(res, ["physical", "stable"], []))
+
+    steadystates = reduce(hcat, stable[index])
+    distances    = vec(sum(abs2.(steadystates .- state), dims=1))
+    return argmin(replace(distances, NaN => Inf))
+end
+
+"""Return the indexes and values following stable branches along a 1D sweep.
+When a no stable solutions are found (e.g. in a bifurcation), the next stable solution is calculated by time evolving the previous solution (quench).
+Keyword arguments
+- `y`:  Dependent variable expression (parsed into Symbolics.jl) to evaluate the followed solution branches on .
+- `sweep`: Direction for the sweeping of solutions. A `right` (`left`) sweep proceeds from the first (last) solution, ordered as the sweeping parameter.
+- `tf`: time to reach steady
+- `ϵ`: small random perturbation applied to quenched solution, in a bifurcation in order to favour convergence in cases where multiple solutions are identically accessible (e.g. symmetry breaking into two equal amplitude states)
+"""
+function follow_branch(starting_branch::Int64, res::Result; y="u1^2+v1^2", sweep="right", tf=10000, ϵ=1e-4)
+    sweep_directions = ["left", "right"]
+    sweep ∈ sweep_directions  || error("Only the following (1D) sweeping directions are allowed:  ", sweep_directions)
+
+    # get stable solutions
+    Y  = transform_solutions(res, y)
+    Ys = HarmonicBalance._apply_mask(Y, HarmonicBalance._get_mask(res, ["physical", "stable"], []))
+    Ys = sweep == "left" ? reverse(Ys) : Ys
+
+    followed_branch    = zeros(Int64,length(Y))  # followed branch indexes
+    followed_branch[1] = starting_branch
+
+    p1 = first(keys(res.swept_parameters)) # parameter values
+
+    for i in 2:length(Ys)
+        s = Ys[i][followed_branch[i-1]] # solution amplitude in the current branch and current parameter index
+        if !isnan(s) # the solution is not unstable or unphysical
+            followed_branch[i] = followed_branch[i-1]
+        else # bifurcation found
+            next_index = sweep == "right" ? i : length(Ys)-i+1
+
+            # create a synthetic starting point out of an unphysical solution: quench and time evolve
+            # the actual solution is complex there, i.e. non physical. Take real part for the quench.
+            sol_dict  = get_single_solution(res, branch=followed_branch[i-1], index=next_index)
+
+            var = res.problem.variables
+            var_values_noise = real.(getindex.(Ref(sol_dict), var)) .+ 0.0im .+ ϵ*rand(length(var))
+            for (i, v) in enumerate(var)
+                sol_dict[v] = var_values_noise[i]
+            end
+
+            problem_t = OrdinaryDiffEq.ODEProblem(res.problem.eom, sol_dict, timespan=(0, tf))
+            res_t     = OrdinaryDiffEq.solve(problem_t, OrdinaryDiffEq.Tsit5(), saveat=tf)
+
+            followed_branch[i] = _closest_branch_index(res, res_t.u[end], next_index) # closest branch to final state
+
+            @info "bifurcation @ $p1 = $(real(sol_dict[p1])): switched branch $(followed_branch[i-1]) ➡ $(followed_branch[i])"
+        end
+    end
+    if sweep == "left"
+        Ys = reverse(Ys)
+        followed_branch = reverse(followed_branch)
+    end
+
+    return followed_branch, Ys
+end
+
+
+"""1D plot with the followed branch highlighted"""
+function plot_1D_solutions_branch(starting_branch::Int64, res::Result;
+    x::String, y::String, sweep="right", tf=10000, ϵ=1e-4, class="default", not_class=[], kwargs...)
+    p = plot(res; x=x, y=y, class=class, not_class=not_class, kwargs...)
+
+    followed_branch, Ys = follow_branch(starting_branch, res, y=y, sweep=sweep, tf=tf, ϵ=ϵ)
+    Y_followed = [Ys[param_idx][branch] for (param_idx,branch) in enumerate(followed_branch)]
+    X = real.(res.swept_parameters[HarmonicBalance._parse_expression(x)])
+
+    Plots.plot!(p, X, real.(Y_followed); linestyle=:dash, c=:gray, label = sweep*" sweep", HarmonicBalance._set_Plots_default..., kwargs...)
+    return p
+end
diff --git a/src/plotting_Plots.jl b/src/plotting_Plots.jl
index 7bd8fabd..bfe74671 100644
--- a/src/plotting_Plots.jl
+++ b/src/plotting_Plots.jl
@@ -94,6 +94,14 @@ function _apply_mask(solns::Array{Vector{ComplexF64}},  booleans)
     factors = replace.(booleans, 0 => NaN)
     map(.*, solns, factors)
 end
+function _apply_mask(solns::Vector{Vector{Vector{ComplexF64}}}, booleans)
+    Nan_vector = NaN.*similar(solns[1][1])
+    new_solns = [
+        [booleans[i][j] ? solns[i][j] : Nan_vector for j in eachindex(solns[i])]
+        for i in eachindex(solns)
+    ]
+    return new_solns
+end
 
 
 """ Project the array `a` into the real axis, warning if its contents are complex. """
diff --git a/src/types.jl b/src/types.jl
index 2b550ac0..99ed0672 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -112,6 +112,7 @@ mutable struct HarmonicEquation
 
     # use a self-referential constructor with _parameters
     HarmonicEquation(equations, variables, nat_eq) = (x = new(equations, variables, Vector{Num}([]), nat_eq); x.parameters=_parameters(x); x)
+    HarmonicEquation(equations, variables, parameters, natural_equation) = new(equations, variables, parameters, natural_equation)
 end
 
 
diff --git a/test/hysteresis_sweep.jl b/test/hysteresis_sweep.jl
new file mode 100644
index 00000000..dd1786c8
--- /dev/null
+++ b/test/hysteresis_sweep.jl
@@ -0,0 +1,24 @@
+using HarmonicBalance, OrdinaryDiffEq
+
+@variables α ω ω0 F γ η t x(t); # declare constant variables and a function x(t)
+
+diff_eq = DifferentialEquation(d(x,t,2) + ω0*x + α*x^3 + γ*d(x,t) + η*x^2*d(x,t) ~ F*cos(ω*t), x) # define ODE
+
+add_harmonic!(diff_eq, x, ω) # specify the ansatz x = u(T) cos(ωt) + v(T) sin(ωt)
+harmonic_eq = get_harmonic_equations(diff_eq) # implement ansatz to get harmonic equations
+
+fixed = (α => 1, ω0 => 1.0, γ => 0.005, F => 0.005, η => 0.2)   # fixed parameters
+varied = ω => range(0.95, 1.1, 10)           # range of parameter values
+result = get_steady_states(harmonic_eq, varied, fixed, show_progress=false)
+
+followed_branches, _ = follow_branch(1, result)
+
+@test first(followed_branches) ≠ last(followed_branches)
+
+plot_1D_solutions_branch(1, result, x="ω", y="√(u1^2+v1^2)")
+
+plot_linear_response(result, x, followed_branches, Ω_range=range(0.9,1.1,10), show=false);
+
+# Check if an empty brspectrum can be plotted
+followed_branches[6] = 3
+plot_linear_response(result, x, followed_branches, Ω_range=range(0.9,1.1,10), show=false, force=true);
diff --git a/test/linear_response.jl b/test/linear_response.jl
index 001f2ce9..3b427137 100644
--- a/test/linear_response.jl
+++ b/test/linear_response.jl
@@ -1,4 +1,5 @@
 using HarmonicBalance
+
 import HarmonicBalance.LinearResponse.plot_linear_response
 import HarmonicBalance.LinearResponse.plot_rotframe_jacobian_response
 import HarmonicBalance.LinearResponse.plot_eigenvalues
@@ -11,8 +12,10 @@ harmonic_eq = get_harmonic_equations(diff_eq)
 
 fixed = (α => 1, ω0 => 1.0, γ => 1e-2, F => 1e-6)
 varied = ω => range(0.9, 1.1, 10)
-result = get_steady_states(harmonic_eq, varied, fixed)
 
-plot_linear_response(result, x, branch=1, Ω_range=range(0.9,1.1,10), order=1, logscale=true);
-plot_rotframe_jacobian_response(result, Ω_range=range(-0.2,0.2,10), branch=1, logscale=true);
-plot_eigenvalues(result, branch=1)
+result = get_steady_states(harmonic_eq, varied, fixed, show_progress=false)
+
+plot_linear_response(result, x, branch=1, Ω_range=range(0.9,1.1,10), order=1, logscale=true)
+
+plot_rotframe_jacobian_response(result, Ω_range=range(0.01,1.1,10), branch=1, logscale=true)
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 9c4e8969..3b476427 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,7 +13,8 @@ files = [
     "transform_solutions.jl",
     "plotting.jl",
     "time_evolution.jl",
-    "krylov.jl"
+    "krylov.jl",
+    "hysteresis_sweep.jl"
     ]
 
 for file in files