diff --git a/Project.toml b/Project.toml
index 71f9b6375..bd71154a5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Symbolics"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 authors = ["Shashi Gowda <gowda@mit.edu>"]
-version = "6.6.1"
+version = "6.7.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/docs/src/manual/solver.md b/docs/src/manual/solver.md
index 7b4a33013..3380d6632 100644
--- a/docs/src/manual/solver.md
+++ b/docs/src/manual/solver.md
@@ -63,6 +63,14 @@ Symbolics.symbolic_solve(eqs, [x,y,z])
 - [ ] Systems of polynomial equations with parameters and positive dimensional systems
 - [ ] Inequalities
 
+### Expressions we can not solve (but aim to)
+```
+# Mathematica
+
+In[1]:= Reduce[x^2 - x - 6 > 0, x]
+Out[1]= x < -2 || x > 3
+```
+
 # References
 
 [^1]: [Rouillier, F. Solving Zero-Dimensional Systems Through the Rational Univariate Representation. AAECC 9, 433–461 (1999).](https://doi.org/10.1007/s002000050114)
diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl
index c8648f3f4..ebf5174a5 100644
--- a/ext/SymbolicsGroebnerExt.jl
+++ b/ext/SymbolicsGroebnerExt.jl
@@ -108,6 +108,7 @@ end
 # Given a GB in k[params][vars] produces a GB in k(params)[vars]
 function demote(gb, vars::Vector{Num}, params::Vector{Num})
     isequal(gb, [1]) && return gb 
+
     gb = Symbolics.wrap.(SymbolicUtils.toterm.(gb))
     Symbolics.check_polynomial.(gb)
 
@@ -126,7 +127,7 @@ function demote(gb, vars::Vector{Num}, params::Vector{Num})
     ring_param, params_demoted = Nemo.polynomial_ring(Nemo.base_ring(ring_flat), map(string, nemo_params))
     ring_demoted, vars_demoted = Nemo.polynomial_ring(Nemo.fraction_field(ring_param), map(string, nemo_vars), internal_ordering=:lex)
     varmap = Dict((nemo_vars .=> vars_demoted)..., (nemo_params .=> params_demoted)...)
-    gb_demoted = map(f -> nemo_crude_evaluate(f, varmap), nemo_gb)
+    gb_demoted = map(f -> ring_demoted(nemo_crude_evaluate(f, varmap)), nemo_gb)
     result = empty(gb_demoted)
     while true
         gb_demoted = map(f -> Nemo.map_coefficients(c -> c // Nemo.leading_coefficient(f), f), gb_demoted)
@@ -176,6 +177,7 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
     # Use a new variable to separate the input polynomials (Reference above)
     new_var = gen_separating_var(vars)
     old_len = length(vars)
+    old_vars = deepcopy(vars)
     vars = vcat(vars, new_var)
 
     new_eqs = []
@@ -204,6 +206,13 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
             return []
         end
 
+        for i in reverse(eachindex(new_eqs))
+            all_present = Symbolics.get_variables(new_eqs[i])
+            if length(intersect(all_present, vars)) < 1
+                deleteat!(new_eqs, i)
+            end
+        end
+
         new_eqs = demote(new_eqs, vars, params)
         new_eqs = map(Symbolics.unwrap, new_eqs)
 
@@ -233,7 +242,10 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa
         end
 
         # non-cyclic case
-        n_iterations > 10 && return []
+        if n_iterations > 10 
+            warns && @warn("symbolic_solve can not currently solve this system of polynomials.")
+            return nothing
+        end
 
         n_iterations += 1
     end
@@ -295,11 +307,13 @@ function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplici
     isempty(tr_basis) && return nothing
     vars_gen = setdiff(vars, tr_basis)
     sol = solve_zerodim(eqs, vars_gen; dropmultiplicity=dropmultiplicity, warns=warns)
+
     for roots in sol
         for x in tr_basis
             roots[x] = x
         end
     end
+
     sol
 end
 
@@ -313,9 +327,8 @@ PrecompileTools.@setup_workload begin
     PrecompileTools.@compile_workload begin
         symbolic_solve(equation1, x)
         symbolic_solve(equation_actually_polynomial)
-        symbolic_solve(simple_linear_equations, [x, y])
-        symbolic_solve(equations_intersect_sphere_line, [x, y, z])
-        symbolic_solve([x^2 - a^2, x + a], x)
+        symbolic_solve(simple_linear_equations, [x, y], warns=false)
+        symbolic_solve(equations_intersect_sphere_line, [x, y, z], warns=false)
     end
 end
 
diff --git a/src/solver/main.jl b/src/solver/main.jl
index 0704abe72..c2540d77e 100644
--- a/src/solver/main.jl
+++ b/src/solver/main.jl
@@ -247,7 +247,7 @@ implemented in the function `get_roots` and its children.
 # Examples
 
 """
-function solve_univar(expression, x; dropmultiplicity = true)
+function solve_univar(expression, x; dropmultiplicity=true)
     args = []
     mult_n = 1
     expression = unwrap(expression)
diff --git a/test/solver.jl b/test/solver.jl
index fb52760d9..cd2e8c42c 100644
--- a/test/solver.jl
+++ b/test/solver.jl
@@ -236,13 +236,13 @@ end
     # cyclic 3
     @variables z1 z2 z3
     eqs = [z1 + z2 + z3, z1*z2 + z1*z3 + z2*z3, z1*z2*z3 - 1]
-    sol = Symbolics.symbolic_solve(eqs, [z1,z2,z3])
+    sol = symbolic_solve(eqs, [z1,z2,z3])
     backward = [Symbolics.substitute(eqs, s) for s in sol]
     @test all(x -> all(isapprox.(eval(Symbolics.toexpr(x)), 0; atol=1e-6)), backward)
 
     @variables x y
     eqs = [2332//232*x + 2131232*y - 1//343434, x + y + 1]
-    sol = Symbolics.symbolic_solve(eqs, [x,y])
+    sol = symbolic_solve(eqs, [x,y])
     backward = [Symbolics.substitute(eqs, s) for s in sol]
     @test all(x -> all(isapprox.(eval(Symbolics.toexpr(x)), 0; atol=1e-6)), backward)
 
@@ -259,10 +259,12 @@ end
         # at most 4 roots by Bézout's theorem
         rand_eq(xs, d) = rand(-10:10) + rand(-10:10)*x + rand(-10:10)*y + rand(-10:10)*x*y + rand(-10:10)*x^2 + rand(-10:10)*y^2
         eqs = [rand_eq([x,y],2), rand_eq([x,y],2)]
-        sol = Symbolics.symbolic_solve(eqs, [x,y])
+        sol = symbolic_solve(eqs, [x,y])
         backward = [Symbolics.substitute(eqs, s) for s in sol]
         @test all(x -> all(isapprox.(eval(Symbolics.toexpr(x)), 0; atol=1e-6)), backward)
     end
+
+    @test isnothing(symbolic_solve([x^2, x*y, y^2], [x,y], warns=false))
 end
 
 @testset "Multivar parametric" begin
@@ -277,8 +279,17 @@ end
     @test isnothing(symbolic_solve([x*y - a, sin(x)], [x, y]))
 
     @variables t w u v
-    sol = symbolic_solve([t*w - 1 ~ 4, u + v + w ~ 1], [t,w,u,v])
-    @test isequal(sol, [Dict(u => u, t => -5 / (-1 + u + v), v => v, w => 1 - u - v)])
+    sol = symbolic_solve([t*w - 1 ~ 4, u + v + w ~ 1], [t,w])
+    @test isequal(sol, [Dict(t => -5 / (-1 + u + v), w => 1 - u - v)])
+
+    sol = symbolic_solve([x-y, y-z], [x])
+    @test isequal(sol, [Dict(x=>z)])
+
+    sol = symbolic_solve([x-y, y-z], [x, y])
+    @test isequal(sol, [Dict(x=>z, y=>z)])
+
+    sol = symbolic_solve([x + y - z, y - z], [x])
+    @test isequal(sol, [Dict(x=>0)])
 end
 
 @testset "Factorisation" begin