diff --git a/src/LinearResponse/plotting.jl b/src/LinearResponse/plotting.jl
index 444703cd..a947aa90 100644
--- a/src/LinearResponse/plotting.jl
+++ b/src/LinearResponse/plotting.jl
@@ -138,7 +138,7 @@ function plot_linear_response(
     kwargs...,
 )
     length(size(res.solutions)) != 1 &&
-        error("1D plots of not-1D datasets are usually a bad idea.")
+        error("The results are two dimensional. Consider using the `cut` keyword.")
     stable = classify_branch(res, branch, "stable") # boolean array
 
     X = Vector{Float64}(collect(values(res.swept_parameters))[1][stable])
@@ -177,7 +177,7 @@ function plot_linear_response(
     kwargs...,
 )
     length(size(res.solutions)) != 1 &&
-        error("1D plots of not-1D datasets are usually a bad idea.")
+        error("The results are two dimensional. Consider using the `cut` keyword.")
 
     X = Vector{Float64}(collect(first(values(res.swept_parameters))))
 
@@ -230,7 +230,7 @@ function plot_rotframe_jacobian_response(
     kwargs...,
 )
     length(size(res.solutions)) != 1 &&
-        error("1D plots of not-1D datasets are usually a bad idea.")
+        error("The results are two dimensional. Consider using the `cut` keyword.")
     stable = classify_branch(res, branch, "stable") # boolean array
 
     Ω_range = vcat(Ω_range)
@@ -277,14 +277,22 @@ function plot_eigenvalues(
     filter = _get_mask(res, class)
     filter_branch = map(x -> getindex(x, branch), replace.(filter, 0 => NaN))
 
-    dim(res) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
+    dim(res) != 1 &&
+        error("The results are two dimensional. Consider using the `cut` keyword.")
     x = string(first(keys(res.swept_parameters)))
     varied = Vector{Float64}(collect(first(values(res.swept_parameters))))
 
-    eigenvalues = [
-        eigvals(res.jacobian(get_single_solution(res; branch=branch, index=i))) for
-        i in eachindex(varied)
-    ]
+    eigenvalues = map(eachindex(varied)) do i
+        jac = res.jacobian(get_single_solution(res; branch=branch, index=i))
+        if any(isnan, jac)
+            throw(
+                ErrorException(
+                    "The branch contains NaN values. Likely, the branch has non-physical solutions in the parameter sweep",
+                ),
+            )
+        end
+        eigvals(jac)
+    end
     eigenvalues_filtered = map(.*, eigenvalues, filter_branch)
 
     eigenvectors = [
diff --git a/src/plotting_Plots.jl b/src/plotting_Plots.jl
index b7432940..572b6f1c 100644
--- a/src/plotting_Plots.jl
+++ b/src/plotting_Plots.jl
@@ -159,7 +159,8 @@ function plot1D(
         end
     end
 
-    dim(res) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
+    dim(res) != 1 &&
+        error("The results are two dimensional. Consider using the `cut` keyword.")
     x = x == "default" ? string(first(keys(res.swept_parameters))) : x
     X = transform_solutions(res, x; branches=branches)
     Y = transform_solutions(res, y; branches=branches, realify=true)
diff --git a/test/linear_response.jl b/test/linear_response.jl
index 6b961c54..c0d86f55 100644
--- a/test/linear_response.jl
+++ b/test/linear_response.jl
@@ -23,3 +23,21 @@ plot_rotframe_jacobian_response(
 
 plot_eigenvalues(result; branch=1)
 plot_eigenvalues(result; branch=1, type=:re, class="all")
+
+@testset begin
+    @variables α λ ω0 ω ωₚ F t x(t)
+    diff_eq = DifferentialEquation(
+        d(x, t, 2) + (ω0^2 - λ * cos(2 * ω * t)) * x + α * x^3 + γ * d(x, t) ~
+            F * cos(ωₚ * t),
+        x,
+    )
+
+    add_harmonic!(diff_eq, x, ω)
+    add_harmonic!(diff_eq, x, ωₚ)
+    harmonic_eq = get_harmonic_equations(diff_eq)
+
+    fixed = (γ => 0.008, ω0 => 1.0, α => 1.0, F => 0.0, ωₚ => 1.0, λ => 0.016)
+    varied = (ω => range(0.995, 1.005, 20))
+    result_ω = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED)
+    @test_throws ErrorException plot_eigenvalues(result_ω; branch=1)
+end