diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl
index dd767c2b6..a7bd4fcaa 100644
--- a/src/RationalFunctionFields/IdealMQS.jl
+++ b/src/RationalFunctionFields/IdealMQS.jl
@@ -194,9 +194,8 @@ function fractionfree_generators_raw(mqs::IdealMQS)
     K = base_ring(ring_params)
     varnames = map(string, Nemo.symbols(ring_params))
     # The hope is that new variables' names would not intersect with the old ones
-    @assert mqs.sat_var_index == length(varnames) + 1
     old_varnames = map(i -> "y$i", 1:length(varnames))
-    new_varnames = map(i -> "라$i", 1:(length(varnames) + 1))
+    new_varnames = map(i -> "__var_$i", 1:(length(varnames) + 1))
     if !isempty(intersect(old_varnames, new_varnames))
         @warn "Intersection in two sets of variables! $varnames $new_varnames"
     end
@@ -221,7 +220,7 @@ function fractionfree_generators_raw(mqs::IdealMQS)
     end
     polys[end] = sat_y
     main_var_indices = 1:(length(varnames) + 1)
-    param_var_indices = (main_var_indices + 1):length(big_vars)
+    param_var_indices = (length(varnames) + 2):length(big_vars)
     return polys, main_var_indices, param_var_indices
 end
 
diff --git a/test/RationalFunctionFields/raw_generators.jl b/test/RationalFunctionFields/raw_generators.jl
new file mode 100644
index 000000000..e41c114b3
--- /dev/null
+++ b/test/RationalFunctionFields/raw_generators.jl
@@ -0,0 +1,61 @@
+# Test that the following pipeline works:
+# ODE -> IO equations -> raw identifiable functions -> RFF -> MQS -> (maybe specialize) -> ideal generators
+import Groebner
+
+@testset "Raw Generators of RFF" begin
+    # SIWR
+    ode = StructuralIdentifiability.@ODEmodel(
+        S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
+        I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
+        W'(t) = xi * (I(t) - W(t)),
+        R'(t) = gam * I(t) - (mu + a) * R(t),
+        y(t) = k * I(t)
+    )
+
+    io_eqs = StructuralIdentifiability.find_ioequations(ode)
+    id_funcs, bring = StructuralIdentifiability.extract_identifiable_functions_raw(
+        io_eqs,
+        ode,
+        empty(ode.parameters),
+        true,
+    )
+
+    param_ring, _ = polynomial_ring(
+        base_ring(bring),
+        map(string, ode.parameters),
+        internal_ordering = Nemo.internal_ordering(bring),
+    )
+
+    id_funcs_no_states = map(
+        polys -> map(
+            poly -> StructuralIdentifiability.parent_ring_change(poly, param_ring),
+            polys,
+        ),
+        id_funcs[:no_states],
+    )
+
+    rff = StructuralIdentifiability.RationalFunctionField(id_funcs_no_states)
+
+    # Part 1: mod p and specialized
+    p = Nemo.Native.GF(2^62 + 135)
+    StructuralIdentifiability.ParamPunPam.reduce_mod_p!(rff.mqs, p)
+    point = rand(
+        p,
+        length(Nemo.gens(StructuralIdentifiability.ParamPunPam.parent_params(rff.mqs))),
+    )
+    eqs = StructuralIdentifiability.ParamPunPam.specialize_mod_p(rff.mqs, point)
+    gb = Groebner.groebner(eqs, ordering = Groebner.DegRevLex())
+    # GB is linear
+    @test length(gb) == length(gens(parent(eqs[1])))
+    expected = 10e6
+    str = join(map(string, eqs), ",")
+    @info "" length(str)
+    @test abs(length(str) - expected) / expected * 100 < 5
+
+    # Part 2: over Q
+    eqs = StructuralIdentifiability.fractionfree_generators_raw(rff.mqs)[1]
+    expected = 26e6
+    str = join(map(string, eqs), ",")
+    @info "" length(str)
+    @test abs(length(str) - expected) / expected * 100 < 5
+end