diff --git a/src/algorithms/dev.ipynb b/src/algorithms/dev.ipynb
index c8f286d..90c9ffb 100644
--- a/src/algorithms/dev.ipynb
+++ b/src/algorithms/dev.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 75,
+   "execution_count": 146,
    "metadata": {},
    "outputs": [
     {
@@ -23,7 +23,7 @@
     "    R, variables = AlgebraicSolving.polynomial_ring(\n",
     "        AlgebraicSolving.GF(101),[\"x\",\"y\",\"u\",\"v\",\"l\",\"dl\"], \n",
     "        internal_ordering=:degrevlex)\n",
-    "    (x,y,u,v,l,dl) = variables\n",
+    "    (x,y,u,v,l,dl) = R_vars\n",
     "\n",
     "    derivatives = Dict(\n",
     "        x => u,\n",
@@ -33,19 +33,19 @@
     "        l => dl)\n",
     "\n",
     "    ideal = AlgebraicSolving.Ideal([x^2 + y^2 - 1])\n",
-    "    return ideal, derivatives, R\n",
+    "    return ideal, derivatives, R, R_vars\n",
     "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 80,
+   "execution_count": 154,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "_manage_rings (generic function with 1 method)"
+       "_swap_vars"
       ]
      },
      "metadata": {},
@@ -53,63 +53,145 @@
     }
    ],
    "source": [
-    "function _manage_rings(ideal, derivatives, R)\n",
+    "using AbstractAlgebra\n",
+    "using AlgebraicSolving\n",
+    "\n",
+    "\"\"\"\n",
+    "    diff_op(q, derivatives)\n",
+    "\n",
+    "    This function evaluates the linear differential operator from \n",
+    "    Definition 7 and Notation 3 in the thesis for a polynomial and derivatives.\n",
+    "    The inputs should not explicitly depend on time.\n",
+    "\n",
+    "    # Arguments\n",
+    "    - `q`: a polynomial\n",
+    "    - `derivatives`: a dictionary of derivatives\n",
+    "\n",
+    "    # Returns\n",
+    "    - the linear differential operator applied to the polynomial\n",
+    "\"\"\"\n",
+    "function diff_op(q, derivatives)    \n",
+    "    n = length(q.parent.data.S)\n",
+    "    @assert length(derivatives) <= n \"There can't be more derivatives than variables.\"\n",
+    "    \n",
+    "    result = 0\n",
+    "    for (var, value) in derivatives\n",
+    "        result += value * AbstractAlgebra.derivative(q, var)\n",
+    "    end\n",
+    "    return result\n",
+    "end\n",
+    "\n",
+    "\"\"\"\n",
+    "    intersect(G, S_vars)\n",
+    "\n",
+    "    This function computes the intersection of a Groebner basis \n",
+    "    with a set of variables. It is important that a proper elimination\n",
+    "    order is used to ensure that the variables that are to be eliminated\n",
+    "    by intersection can actually be eliminated.\n",
+    "\n",
+    "    An example would be lexico-graphical ordering where the variables that are\n",
+    "    to be eliminated are the first ones. Or alternatively and superior in\n",
+    "    practice is to use block ordering where the variables that are to be\n",
+    "    eliminated are in the larger block w.r.t. to the ordering.\n",
+    "\n",
+    "    # Arguments\n",
+    "    - `G`: a Groebner basis w.r.t a monomial odering that eliminates all \n",
+    "       variables except the ones in `S_vars`\n",
+    "    - `S_vars`: the list of variables that are not to be eliminated\n",
+    "\n",
+    "    # Returns\n",
+    "    - the intersection of the Groebner basis with the subring S\n",
+    "\"\"\"\n",
+    "function intersect(G, S_vars)\n",
+    "    sub_ideal = []\n",
+    "    for generator in G\n",
+    "        symbols = [Symbol(var) for var in AbstractAlgebra.vars(generator)]\n",
+    "        if issubset(symbols, S_vars)\n",
+    "            push!(sub_ideal, generator)\n",
+    "        end\n",
+    "    end\n",
+    "    return sub_ideal\n",
+    "end\n",
+    "\n",
+    "\"\"\"\n",
+    "    _manage_rings(ideal, derivatives, R)\n",
+    "\n",
+    "    This function creates a new ring that permits the elimination of variables.\n",
+    "\n",
+    "    # Arguments\n",
+    "    - `ideal`: an ideal\n",
+    "    - `derivatives`: a dictionary of derivatives\n",
+    "    - `R`: a polynomial ring with two blocks of variables, where the first \n",
+    "       block corresponds to the variables in `R` that are not in `derivatives` \n",
+    "       and the second block corresponds to the variables in `derivatives`\n",
+    "\n",
+    "    # Returns\n",
+    "    - the new ring\n",
+    "    - the new variables\n",
+    "    - a map from the old variables to the new variables\n",
+    "\"\"\"\n",
+    "function _manage_rings(derivatives, R)\n",
     "    S_proper = []\n",
     "    R_elim = []\n",
+    "    s = length(R.data.S) - length(derivatives)\n",
+    "    map = Dict()\n",
+    "    \n",
+    "    index_n_el = s + 1\n",
+    "    index_el = 1\n",
+    "\n",
+    "    # Move the variables to the correct blocks and create a map\n",
     "    for var in R.data.S\n",
     "        if Symbol(var) in [Symbol(key) for key in keys(derivatives)]\n",
     "            push!(S_proper, Symbol(var))\n",
+    "            map[Symbol(var)] = index_n_el\n",
+    "            index_n_el += 1\n",
     "        else    \n",
     "            push!(R_elim, Symbol(var))\n",
+    "            map[Symbol(var)] = index_el\n",
+    "            index_el += 1\n",
     "        end\n",
     "    end\n",
     "\n",
     "    # Create new ring, where the first block is the eliminated variables\n",
     "    R_vars_proper = append!(R_elim, S_proper)\n",
-    "    R_names = [string(var) for var in R_vars_proper]\n",
-    "    R_new, R_names = AlgebraicSolving.polynomial_ring(\n",
-    "        base_ring(R), R_names, internal_ordering=:degrevlex)\n",
+    "    R_new_vars = [string(var) for var in R_vars_proper]\n",
+    "    R_new, R_new_vars = AlgebraicSolving.polynomial_ring(\n",
+    "        base_ring(R), R_new_vars, internal_ordering=:degrevlex)\n",
     "    \n",
-    "    return R_new\n",
+    "    return R_new, R_new_vars, map\n",
+    "end\n",
+    "\n",
+    "\n",
+    "\"\"\"\n",
+    "    _swap_vars(poly, R_new, R_new_vars, map)\n",
+    "\n",
+    "    This function swaps the variables in a polynomial to a new ring.\n",
+    "\n",
+    "    # Arguments\n",
+    "    - `poly`: a polynomial\n",
+    "    - `R_new`: a new ring\n",
+    "    - `R_new_vars`: the new variables\n",
+    "    - `map`: a map from the old variables to the new variables\n",
+    "\n",
+    "    # Returns\n",
+    "    - the polynomial in the new ring\n",
+    "\"\"\"\n",
+    "function _swap_vars(poly, R_vars, R_new_vars, map)\n",
+    "    vars_subst = [R_new_vars[map[Symbol(var)]] for var in R_vars]\n",
+    "    poly_new = poly(vars_subst...)\n",
+    "    return poly_new\n",
     "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 86,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Multivariate polynomial ring in 6 variables over GF(101)[:x, :y, :u, :v, :l, :dl]\n",
-      "Multivariate polynomial ring in 6 variables over GF(101)[:dl, :x, :y, :u, :v, :l]\n"
-     ]
-    }
-   ],
-   "source": [
-    "ideal, derivatives, R = simple_pendulum()\n",
-    "println(R, R.data.S)\n",
-    "R_new = _manage_rings(ideal, derivatives, R)\n",
-    "println(R_new, R_new.data.S)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 96,
+   "execution_count": 184,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "6-element Vector{Int64}:\n",
-       " 2\n",
-       " 0\n",
-       " 0\n",
-       " 0\n",
-       " 0\n",
-       " 0"
+       "differential_basis"
       ]
      },
      "metadata": {},
@@ -117,25 +199,102 @@
     }
    ],
    "source": [
-    "s = ideal[1]\n",
-    "exponen = AbstractAlgebra.exponent_vector(s, 1)"
+    "\n",
+    "\n",
+    "\"\"\"\n",
+    "    differential_basis(ideal, derivatives, R, nf=false, info_level=0)\n",
+    "\n",
+    "    This function computes the differential Gröbner basis of an ideal or \n",
+    "    an approximation of a full differential ideal if not all derivatives\n",
+    "    are known. See Algorithm 4 in the thesis for reference.\n",
+    "\n",
+    "    # Arguments\n",
+    "    - `ideal`: an ideal\n",
+    "    - `derivatives`: a dictionary of derivatives\n",
+    "    - `R`: a polynomial ring with two blocks of variables, where the first \n",
+    "       block corresponds to the variables in `R` that are not in `derivatives` \n",
+    "       and the second block corresponds to the variables in `derivatives`\n",
+    "    - `R_vars`: the variables in the ring, only for elimination, else use `[]``\n",
+    "    - `nf`: a boolean indicating whether to compute the normal form\n",
+    "    - `info_level`: an integer indicating the level of information to print\n",
+    "\n",
+    "    # Returns\n",
+    "    - the differential Gröbner basis\n",
+    "\"\"\"\n",
+    "function differential_basis(ideal, derivatives, R, R_vars, nf=false, info_level=0)\n",
+    "    n = R.data.nvars\n",
+    "    s = length(derivatives)    \n",
+    "    eliminate = n - s\n",
+    "\n",
+    "    # If elimination is necessary reorganize the ring and substitute the variables\n",
+    "    if eliminate > 0\n",
+    "        R_new, R_new_vars, map = _manage_rings(derivatives, R)\n",
+    "        \n",
+    "        ideal_new_gens = [_swap_vars(elem, R_vars, R_new_vars, map) for elem in ideal.gens]\n",
+    "        ideal = AlgebraicSolving.Ideal(ideal_new_gens)\n",
+    "\n",
+    "        derivatives_new = Dict()\n",
+    "        for (var, expr) in derivatives\n",
+    "            derivatives_new[_swap_vars(var, R_vars, R_new_vars, map)] = _swap_vars(expr, R_vars, R_new_vars, map)\n",
+    "        end\n",
+    "        derivatives = derivatives_new\n",
+    "    end \n",
+    "\n",
+    "    # Infer and create the subring\n",
+    "    S_vars = [Symbol(var) for var in R_new_vars[eliminate+1:end]]\n",
+    "    \n",
+    "    # Start computing the differential basis\n",
+    "    G1 = groebner_basis(ideal)\n",
+    "    pG1 = [diff_op(g, derivatives) for g in intersect(G1, S_vars)]\n",
+    "    if nf\n",
+    "        pG1 = [AlgebraicSolving.normal_form(pg, AlgebraicSolving.Ideal(G1)) for pg in pG1]\n",
+    "    end\n",
+    "    append!(pG1, G1)\n",
+    "    G2 = groebner_basis(AlgebraicSolving.Ideal(pG1), eliminate=eliminate,\n",
+    "                        intersect=false, info_level=info_level)\n",
+    "    if info_level > 0\n",
+    "        i = 1\n",
+    "        println(\"iteration \", i)\n",
+    "        println(\"#G = \", length(G1))\n",
+    "    end\n",
+    "\n",
+    "    # Repeat until closed under diff_op\n",
+    "    while G1 != G2\n",
+    "        if info_level > 0\n",
+    "            i += 1\n",
+    "            println(\"iteration \", i)\n",
+    "            println(\"#G = \", length(G2))\n",
+    "        end\n",
+    "        G1 = G2\n",
+    "        pG1 = [diff_op(g, derivatives) for g in intersect(G1, S_vars)]\n",
+    "        if nf == true\n",
+    "            pG1 = [AlgebraicSolving.normal_form(pg, AlgebraicSolving.Ideal(G1)) for pg in pG1]\n",
+    "        end\n",
+    "            append!(pG1, G1)\n",
+    "        G2 = groebner_basis(AlgebraicSolving.Ideal(pG1), eliminate=eliminate,\n",
+    "                            intersect=false, info_level=info_level)\n",
+    "    end\n",
+    "    return G1\n",
+    "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 97,
+   "execution_count": 186,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "6-element Vector{Symbol}:\n",
-       " :dl\n",
-       " :x\n",
-       " :y\n",
-       " :u\n",
-       " :v\n",
-       " :l"
+       "8-element Vector{FqMPolyRingElem}:\n",
+       " u^2 + v^2 + 100*y + l\n",
+       " x*u + y*v\n",
+       " x^2 + y^2 + 100\n",
+       " y*u*v + 100*x*v^2 + x*y + 100*x*l\n",
+       " y^2*u + 100*x*y*v + 100*u\n",
+       " y^3 + 100*y^2*l + v^2 + 100*y + l\n",
+       " x*y^2 + 100*x*y*l + u*v\n",
+       " dl + 98*v"
       ]
      },
      "metadata": {},
@@ -143,7 +302,8 @@
     }
    ],
    "source": [
-    "R_new.data.S"
+    "ideal, derivatives, R, R_vars = simple_pendulum()\n",
+    "res = differential_basis(ideal, derivatives, R, R_vars)"
    ]
   },
   {
diff --git a/src/algorithms/main.jl b/src/algorithms/main.jl
index 37e0621..19e0434 100644
--- a/src/algorithms/main.jl
+++ b/src/algorithms/main.jl
@@ -58,6 +58,75 @@ function intersect(G, S_vars)
     return sub_ideal
 end
 
+"""
+    _manage_rings(ideal, derivatives, R)
+
+    This function creates a new ring that permits the elimination of variables.
+
+    # Arguments
+    - `ideal`: an ideal
+    - `derivatives`: a dictionary of derivatives
+    - `R`: a polynomial ring with two blocks of variables, where the first 
+       block corresponds to the variables in `R` that are not in `derivatives` 
+       and the second block corresponds to the variables in `derivatives`
+
+    # Returns
+    - the new ring
+    - the new variables
+    - a map from the old variables to the new variables
+"""
+function _manage_rings(derivatives, R)
+    S_proper = []
+    R_elim = []
+    s = length(R.data.S) - length(derivatives)
+    map = Dict()
+    
+    index_n_el = s + 1
+    index_el = 1
+
+    # Move the variables to the correct blocks and create a map
+    for var in R.data.S
+        if Symbol(var) in [Symbol(key) for key in keys(derivatives)]
+            push!(S_proper, Symbol(var))
+            map[Symbol(var)] = index_n_el
+            index_n_el += 1
+        else    
+            push!(R_elim, Symbol(var))
+            map[Symbol(var)] = index_el
+            index_el += 1
+        end
+    end
+
+    # Create new ring, where the first block is the eliminated variables
+    R_vars_proper = append!(R_elim, S_proper)
+    R_new_vars = [string(var) for var in R_vars_proper]
+    R_new, R_new_vars = AlgebraicSolving.polynomial_ring(
+        base_ring(R), R_new_vars, internal_ordering=:degrevlex)
+    
+    return R_new, R_new_vars, map
+end
+
+
+"""
+    _swap_vars(poly, R_new, R_new_vars, map)
+
+    This function swaps the variables in a polynomial to a new ring.
+
+    # Arguments
+    - `poly`: a polynomial
+    - `R_new`: a new ring
+    - `R_new_vars`: the new variables
+    - `map`: a map from the old variables to the new variables
+
+    # Returns
+    - the polynomial in the new ring
+"""
+function _swap_vars(poly, R_vars, R_new_vars, map)
+    vars_subst = [R_new_vars[map[Symbol(var)]] for var in R_vars]
+    poly_new = poly(vars_subst...)
+    return poly_new
+end
+
 
 """
     differential_basis(ideal, derivatives, R, nf=false, info_level=0)
@@ -72,18 +141,34 @@ end
     - `R`: a polynomial ring with two blocks of variables, where the first 
        block corresponds to the variables in `R` that are not in `derivatives` 
        and the second block corresponds to the variables in `derivatives`
+    - `R_vars`: the variables in the ring, only needed for elimination
     - `nf`: a boolean indicating whether to compute the normal form
     - `info_level`: an integer indicating the level of information to print
 
     # Returns
     - the differential Gröbner basis
 """
-function differential_basis(ideal, derivatives, R, nf=false, info_level=0)
-    # Infer and create the subring
+function differential_basis(ideal, derivatives, R, R_vars = [], nf=false, info_level=0)
     n = R.data.nvars
-    S_vars = [Symbol(var.first) for var in derivatives]
-    k = length(S_vars)
-    eliminate = n - k
+    s = length(derivatives)    
+    eliminate = n - s
+
+    # If elimination is necessary reorganize the ring and substitute the variables
+    if eliminate > 0
+        R_new, R_new_vars, map = _manage_rings(derivatives, R)
+        
+        ideal_new_gens = [_swap_vars(elem, R_vars, R_new_vars, map) for elem in ideal.gens]
+        ideal = AlgebraicSolving.Ideal(ideal_new_gens)
+
+        derivatives_new = Dict()
+        for (var, expr) in derivatives
+            derivatives_new[_swap_vars(var, R_vars, R_new_vars, map)] = _swap_vars(expr, R_vars, R_new_vars, map)
+        end
+        derivatives = derivatives_new
+    end 
+
+    # Infer and create the subring
+    S_vars = [Symbol(var) for var in R_new_vars[eliminate+1:end]]
     
     # Start computing the differential basis
     G1 = groebner_basis(ideal)
diff --git a/src/systems/chemical.jl b/src/systems/chemical.jl
index 49b5243..7deae82 100644
--- a/src/systems/chemical.jl
+++ b/src/systems/chemical.jl
@@ -11,10 +11,11 @@ using AlgebraicSolving
 """
 function chem_1()
     
-    R, (x1, x2, x3, x4, k1, k2, k3, T1, T2) = AlgebraicSolving.polynomial_ring(
+    R, R_vars = AlgebraicSolving.polynomial_ring(
         AlgebraicSolving.GF(101),["x1","x2","x3","x4","k1","k2","k3","T1","T2"], 
         internal_ordering=:degrevlex)
     
+    (x1, x2, x3, x4, k1, k2, k3, T1, T2) = R_vars
     derivatives = Dict(
         x1 => - k1*x1 + k2*x2*x3,    
         x2 => k1*x1 - k2*x2*x3,
@@ -23,7 +24,7 @@ function chem_1()
         k1 => 0, k2 => 0, k3 => 0, T1 => 0, T2 => 0)
 
     ideal = AlgebraicSolving.Ideal([x1 + x2 - T1, x3 + x4 - T2])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
 
 """
@@ -35,7 +36,7 @@ end
     The example can be found as example 2.5.5 in the thesis.
 """
 function akzo_nobel()
-    R, variables = AlgebraicSolving.polynomial_ring(
+    R, R_vars = AlgebraicSolving.polynomial_ring(
         AlgebraicSolving.GF(101),
         ["dx6","dx7",
         "x1","x2","x3","x4","x5","x6","x7",
@@ -45,7 +46,7 @@ function akzo_nobel()
     # 1/H => H, 1/K => K
     (dx6, dx7, 
     x1, x2, x3, x4, x5, x6, x7, 
-    k1, k2, k3, k4, K, klA, Ks, pCO2, H) = variables
+    k1, k2, k3, k4, K, klA, Ks, pCO2, H) = R_vars
     
     r1 = k1 * x1^4 * x7
     r2 = k2 * x3 * x5
@@ -73,7 +74,7 @@ function akzo_nobel()
         H => 0)
 
     ideal = AlgebraicSolving.Ideal([Ks * x1 * x4 - x6, x7^2 - x2])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
 
 """
@@ -85,7 +86,7 @@ end
     The example can be found as example 2.5.6 in the thesis.
 """
 function fast_slow_reaction()
-    R, variables = AlgebraicSolving.polynomial_ring(AlgebraicSolving.GF(101),
+    R, R_vars = AlgebraicSolving.polynomial_ring(AlgebraicSolving.GF(101),
         ["dRA","dRB","dV2","dCAi","dF","dFi",
         "RA","RB","V1","V2","CA","CAi","CB","CC","F","Fi",
         "Keq","kB"],
@@ -94,7 +95,7 @@ function fast_slow_reaction()
     # 1/Keq => Keq
     (dRA, dRB, dV2, dCAi, dF, dFi, 
     RA, RB, V1, V2, CA, CAi, CB, CC, F, Fi, 
-    Keq, kB) = variables
+    Keq, kB) = R_vars
     
     derivatives = Dict(
         RA => dRA,
@@ -111,5 +112,5 @@ function fast_slow_reaction()
         kB => 0)
 
     ideal = AlgebraicSolving.Ideal([CA - CB * Keq, RB - kB * CB, V1*V2-1])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
diff --git a/src/systems/linear_nn.jl b/src/systems/linear_nn.jl
index 999cdc8..66a2426 100644
--- a/src/systems/linear_nn.jl
+++ b/src/systems/linear_nn.jl
@@ -39,7 +39,7 @@ function linear_nn_2(m, n, r)
         push!(C_flat, elem)
     end
     ideal = AlgebraicSolving.Ideal(C_flat)
-    return ideal, derivatives, R
+    return ideal, derivatives, R, []
 end
 
 """
@@ -86,5 +86,5 @@ function linear_nn_3(m, n, r, s)
         push!(D_flat, elem)
     end
     ideal = AlgebraicSolving.Ideal(D_flat)
-    return ideal, derivatives, R
+    return ideal, derivatives, R, []
 end
diff --git a/src/systems/mechanical.jl b/src/systems/mechanical.jl
index b22f938..7e6cdd4 100644
--- a/src/systems/mechanical.jl
+++ b/src/systems/mechanical.jl
@@ -7,10 +7,10 @@ using AlgebraicSolving
     For reference see Example 2.5.1 in the thesis.
 """
 function simple_pendulum()
-    R, variables = AlgebraicSolving.polynomial_ring(
+    R, R_vars = AlgebraicSolving.polynomial_ring(
         AlgebraicSolving.GF(101),["dl","x","y","u","v","l"], 
         internal_ordering=:degrevlex)
-    (dl,x,y,u,v,l) = variables
+    (dl,x,y,u,v,l) = R_vars
 
     derivatives = Dict(
         x => u,
@@ -20,7 +20,7 @@ function simple_pendulum()
         l => dl)
 
     ideal = AlgebraicSolving.Ideal([x^2 + y^2 - 1])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
 
 """
@@ -30,10 +30,10 @@ end
     For reference see Example 2.5.2 in the thesis.
 """
 function double_pendulum()
-    R, variables = AlgebraicSolving.polynomial_ring(AlgebraicSolving.GF(101),
+    R, R_vars = AlgebraicSolving.polynomial_ring(AlgebraicSolving.GF(101),
         ["dl1","dl2","x1","y1","u1","v1","x2","y2","u2","v2","l1","l2"],
         internal_ordering=:degrevlex)
-    (dl1,dl2,x1,y1,u1,v1,x2,y2,u2,v2,l1,l2) = variables
+    (dl1,dl2,x1,y1,u1,v1,x2,y2,u2,v2,l1,l2) = R_vars
 
     derivatives = Dict(
         x1 => u1,
@@ -48,7 +48,7 @@ function double_pendulum()
         l2 => dl2)
 
     ideal = AlgebraicSolving.Ideal([x1^2 + y1^2 - 1, (x2-x1)^2 + (y2-y1)^2 - 1])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
 
 """
@@ -59,12 +59,12 @@ end
     For reference see Example 2.5.2 in the thesis.
 """
 function triple_pendulum()
-    R, variables = AlgebraicSolving.polynomial_ring(
+    R, R_vars = AlgebraicSolving.polynomial_ring(
         AlgebraicSolving.GF(101),
         ["dl1","dl2","dl3","x1","y1","u1","v1","x2","y2",
         "u2","v2","x3","y3","u3","v3","l1","l2","l3"],
         internal_ordering=:degrevlex)
-    (dl1,dl2,dl3,x1,y1,u1,v1,x2,y2,u2,v2,x3,y3,u3,v3,l1,l2,l3) = variables
+    (dl1,dl2,dl3,x1,y1,u1,v1,x2,y2,u2,v2,x3,y3,u3,v3,l1,l2,l3) = R_vars
     
     derivatives = Dict(
         x1 => u1,
@@ -85,7 +85,7 @@ function triple_pendulum()
 
     ideal = AlgebraicSolving.Ideal(
         [x1^2 + y1^2 - 1, (x2-x1)^2 + (y2-y1)^2 - 1, (x3-x2)^2 + (y3-y2)^2 - 1])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
 
 """
@@ -97,10 +97,11 @@ end
     For reference see Example 2.5.3 in the thesis.
 """
 function masspoint_parabola()
-    R, (dl, p1, p2, p3, v1, v2, v3, l) = AlgebraicSolving.polynomial_ring(
+    R, R_vars = AlgebraicSolving.polynomial_ring(
         AlgebraicSolving.GF(101),["dl","p1","p2","p3","v1","v2","v3","l"], 
         internal_ordering=:degrevlex)
     
+        (dl, p1, p2, p3, v1, v2, v3, l) = R_vars
     derivatives = Dict(    
         p1 => v1,
         p2 => v2,
@@ -111,5 +112,5 @@ function masspoint_parabola()
         l => dl)
 
     ideal = AlgebraicSolving.Ideal([p1^2 + p2^2 - p3])
-    return ideal, derivatives, R
+    return ideal, derivatives, R, R_vars
 end
diff --git a/src/systems/poly_nn.jl b/src/systems/poly_nn.jl
index d2abf42..8b10a0e 100644
--- a/src/systems/poly_nn.jl
+++ b/src/systems/poly_nn.jl
@@ -31,5 +31,5 @@ function poly_nn_2(m, n, r, activation=activation_function)
         push!(C_flat, elem)
     end
     ideal = Ideal(C_flat)
-    return ideal, derivatives, R
+    return ideal, derivatives, R, []
 end
diff --git a/test/algorithms/main.jl b/test/algorithms/main.jl
index 68a2fb8..50b4f23 100644
--- a/test/algorithms/main.jl
+++ b/test/algorithms/main.jl
@@ -38,11 +38,13 @@ end
     using DifferentialBases
     using AlgebraicSolving
 
-    R, (dl,l,v,u,y,x) = AlgebraicSolving.polynomial_ring(
+    R, R_vars = AlgebraicSolving.polynomial_ring(
         AlgebraicSolving.GF(101),
         ["dl","l","v","u","y","x"], 
         internal_ordering=:degrevlex)
 
+    (dl, l, v, u, y, x) = R_vars
+
     derivatives = Dict(
         x => u,
         y => v,
@@ -51,6 +53,8 @@ end
         l => dl
     )
 
+    ideal = AlgebraicSolving.Ideal([x^2 + y^2 - 1])
+
     sol = [
         y^2 + x^2 + 100,
         v*y + u*x,
@@ -62,8 +66,42 @@ end
         dl + 98*v
     ]
 
-    ideal = AlgebraicSolving.Ideal([x^2 + y^2 - 1])
-    @test DifferentialBases.differential_basis(ideal, derivatives, R) == sol
     @test DifferentialBases.differential_basis(
-        AlgebraicSolving.Ideal(sol), derivatives, R, true, 1) == sol
+        ideal, derivatives, R, R_vars) == sol
+    @test DifferentialBases.differential_basis(
+        AlgebraicSolving.Ideal(sol), derivatives, R, R_vars, true, 1) == sol
 end
+
+@testset "Algorithms -> Main -> Ring helpers" begin
+    R, R_vars = AlgebraicSolving.polynomial_ring(
+        AlgebraicSolving.GF(101),
+        ["l","v","u","y","x","dl"], 
+        internal_ordering=:degrevlex)
+
+    (l, v, u, y, x, dl) = R_vars
+
+    derivatives = Dict(
+        x => u,
+        y => v,
+        u => x*l,
+        v => y*l - 1,
+        l => dl
+    )
+
+    ideal = AlgebraicSolving.Ideal([x^2 + y^2 - 1])
+
+    R_new, R_new_vars, map = DifferentialBases._manage_rings(derivatives, R)
+
+    R_1, R_1_vars = AlgebraicSolving.polynomial_ring(
+        AlgebraicSolving.GF(101),
+        ["dl","l","v","u","y","x",], 
+        internal_ordering=:degrevlex)
+    
+    
+    @test R_new == R_1
+
+    ideal_new_gens = [DifferentialBases._swap_vars(elem, R_vars, R_new_vars, map) for elem in ideal.gens]
+    ideal = AlgebraicSolving.Ideal(ideal_new_gens)
+
+    @test parent(ideal[1]) == R_new
+end
\ No newline at end of file
diff --git a/test/systems/chemical.jl b/test/systems/chemical.jl
index f36b14e..b54f497 100644
--- a/test/systems/chemical.jl
+++ b/test/systems/chemical.jl
@@ -1,15 +1,15 @@
 @testset "Systems -> chemical" begin
     using DifferentialBases
 
-    ideal, derivatives, R = DifferentialBases.chem_1()
+    ideal, derivatives, R, R_vars = DifferentialBases.chem_1()
     @test length(ideal.gens) == 2
     @test length(derivatives) == 9
 
-    ideal, derivatives, R = DifferentialBases.akzo_nobel()
+    ideal, derivatives, R, R_vars = DifferentialBases.akzo_nobel()
     @test length(ideal.gens) == 2
     @test length(derivatives) == 16
 
-    ideal, derivatives, R = DifferentialBases.fast_slow_reaction()
+    ideal, derivatives, R, R_vars = DifferentialBases.fast_slow_reaction()
     @test length(ideal.gens) == 3
     @test length(derivatives) == 12
 end
diff --git a/test/systems/linear_nn.jl b/test/systems/linear_nn.jl
index 36dd8c4..25808eb 100644
--- a/test/systems/linear_nn.jl
+++ b/test/systems/linear_nn.jl
@@ -1,22 +1,22 @@
 @testset "Systems -> linear_nn" begin
     using DifferentialBases
 
-    ideal, derivatives, R = DifferentialBases.linear_nn_2(2, 2, 2)
+    ideal, derivatives, R, R_vars = DifferentialBases.linear_nn_2(2, 2, 2)
     @test length(ideal.gens) == 2*2
     @test length(derivatives) == 2*2 + 2*2 + 2*2
     @test any(x -> x != 0, derivatives)
     
-    ideal, derivatives, R = DifferentialBases.linear_nn_2(2, 3, 4)
+    ideal, derivatives, R, R_vars = DifferentialBases.linear_nn_2(2, 3, 4)
     @test length(ideal.gens) == 3*3
     @test length(derivatives) == 2*3 + 3*4 + 2 + 4
     @test any(x -> x != 0, derivatives)
 
-    ideal, derivatives, R = DifferentialBases.linear_nn_3(2, 2, 2, 2)
+    ideal, derivatives, R, R_vars = DifferentialBases.linear_nn_3(2, 2, 2, 2)
     @test length(ideal.gens) == 2*2 + 2*2
     @test length(derivatives) == 2*2 + 2*2 + 2*2 + 2*2
     @test any(x -> x != 0, derivatives)
 
-    ideal, derivatives, R = DifferentialBases.linear_nn_3(2, 3, 4, 5)
+    ideal, derivatives, R, R_vars = DifferentialBases.linear_nn_3(2, 3, 4, 5)
     @test length(ideal.gens) == 3*3 + 4*4
     @test length(derivatives) == 2*3 + 3*4 + 4*5 + 2 + 5
     @test any(x -> x != 0, derivatives)
diff --git a/test/systems/mechanical.jl b/test/systems/mechanical.jl
index 18cbc57..dc254b2 100644
--- a/test/systems/mechanical.jl
+++ b/test/systems/mechanical.jl
@@ -2,19 +2,19 @@
 @testset "Systems -> mechanical" begin
     using DifferentialBases
 
-    ideal, derivatives, R = DifferentialBases.simple_pendulum()
+    ideal, derivatives, R, R_vars = DifferentialBases.simple_pendulum()
     @test length(ideal.gens) == 1
     @test length(derivatives) == 5
 
-    ideal, derivatives, R = DifferentialBases.double_pendulum()
+    ideal, derivatives, R, R_vars = DifferentialBases.double_pendulum()
     @test length(ideal.gens) == 2
     @test length(derivatives) == 10
 
-    ideal, derivatives, R = DifferentialBases.triple_pendulum()
+    ideal, derivatives, R, R_vars = DifferentialBases.triple_pendulum()
     @test length(ideal.gens) == 3
     @test length(derivatives) == 15
 
-    ideal, derivatives, R = DifferentialBases.masspoint_parabola()
+    ideal, derivatives, R, R_vars = DifferentialBases.masspoint_parabola()
     @test length(ideal.gens) == 1
     @test length(derivatives) == 7
 end
diff --git a/test/systems/poly_nn.jl b/test/systems/poly_nn.jl
index f948a9c..385c12a 100644
--- a/test/systems/poly_nn.jl
+++ b/test/systems/poly_nn.jl
@@ -1,12 +1,12 @@
 @testset "Systems -> poly_nn" begin
     using DifferentialBases
 
-    ideal, derivatives, R = DifferentialBases.poly_nn_2(2, 2, 2)
+    ideal, derivatives, R, R_vars = DifferentialBases.poly_nn_2(2, 2, 2)
     @test length(ideal.gens) == 2*2
     @test length(derivatives) == 2*2 + 2*2 + 2*2
     @test any(x -> x != 0, derivatives)
     
-    ideal, derivatives, R = DifferentialBases.poly_nn_2(2, 3, 4)
+    ideal, derivatives, R, R_vars = DifferentialBases.poly_nn_2(2, 3, 4)
     @test length(ideal.gens) == 3*3
     @test length(derivatives) == 2*3 + 3*4 + 2 + 4
     @test any(x -> x != 0, derivatives)