diff --git a/src/ExactODEReduction.jl b/src/ExactODEReduction.jl
index 2f707b3..7c01cd0 100644
--- a/src/ExactODEReduction.jl
+++ b/src/ExactODEReduction.jl
@@ -405,6 +405,28 @@ function find_reductions(
     ChainOfReductions(results)
 end
 
+function find_conservation_laws(system::ODE{P}; seed=nothing, loglevel=Logging.Info) where {P}
+    prev_logger = Logging.global_logger(ConsoleLogger(stderr, loglevel))
+
+    isnothing(seed) && (seed = getnewrandomseed())
+    Random.seed!(seed)
+    @info "Global random seed: $seed"
+
+    eqs = equations_extended(system)
+
+    basis, kerdim, ker, relations = constant_linear_relations(eqs)
+    _relations_sym = ker * basis
+    relations_sym = [_relations_sym[i, :] for i in 1:kerdim]
+
+    Logging.global_logger(prev_logger)
+
+    laws = []
+    for i in 1:kerdim
+        push!(laws, (basis=basis, matrix=relations[i], law=relations_sym[i]))
+    end
+    laws
+end
+
 export ODE, @ODEsystem, equations, vars
 export states, parameters, initial_conditions, parameter_values
 export set_parameter_values, to_state, to_parameter
diff --git a/src/Reduction.jl b/src/Reduction.jl
index b69e51a..e897dcd 100644
--- a/src/Reduction.jl
+++ b/src/Reduction.jl
@@ -38,14 +38,21 @@ struct Reduction{A, B, C}
     new_system::ODE{A}
     new_vars::Dict{A, B}
     old_system::ODE{C}
+    lumping_matrix::Union{Nothing, Matrix{Any}}
 
     function Reduction{A, B, C}(new_system::ODE{A}, new_vars::Dict{A, B}, old_system::ODE{C}) where {A, B, C}
-        return new{A, B, C}(new_system, new_vars, old_system)
+        lumping_matrix = nothing
+        if !isempty(new_vars)
+            lumping_vars = [new_vars[x] for x in gens(parent(new_system))]
+            lumping_combs = [sum(collect(coefficients(comb)) .* exponent_vectors(comb)) for comb in lumping_vars]
+            lumping_matrix = mapreduce(permutedims, vcat, lumping_combs)
+        end
+        return new{A, B, C}(new_system, new_vars, old_system, lumping_matrix)
     end
 
     function Reduction{P}(old_ode::ODE{P}, subspace, parameter_strategy=:inhertiance) where {P}
-        (transformation, new_equations) = perform_change_of_variables(equations_extended(old_ode), subspace)
-    
+        (lumping_matrix, transformation, new_equations) = perform_change_of_variables(equations_extended(old_ode), subspace)
+        
         # (!) assumes the correct order in new_equations, i.e.,
         # ∂y[1] ~ new_equations[1],
         # ∂y[2] ~ new_equations[2],
@@ -82,7 +89,7 @@ struct Reduction{A, B, C}
             @warn "Unknown parameter handling strategy $parameter_strategy, using none"
         end
 
-        return new{elem_type(parent(new_ode)), elem_type(parent(first(transformation))), elem_type(parent(old_ode))}(new_ode, new_vars, old_ode)
+        return new{elem_type(parent(new_ode)), elem_type(parent(first(transformation))), elem_type(parent(old_ode))}(new_ode, new_vars, old_ode, lumping_matrix)
     end
 end
 
diff --git a/src/linalg_utils.jl b/src/linalg_utils.jl
index bbb9ba7..ec4c744 100644
--- a/src/linalg_utils.jl
+++ b/src/linalg_utils.jl
@@ -136,3 +136,24 @@ function construct_jacobians(system)
     ring = parent(first(keys(system)))
     construct_jacobians([system[x] for x in gens(ring)])
 end
+
+#------------------------------------------------------------------------------
+
+# Find the left kernel of a Macaulay matrix of `system`.
+function constant_linear_relations(system::AbstractArray{T}) where {T<:Nemo.MPolyElem}
+    domain = base_ring(first(system))
+    poly_ring = parent(first(system))
+    m = length(system)
+    labels = unique(sort(reduce(vcat, map(collect ∘ monomials, system)), rev=true))
+    n = length(labels)
+    S = Nemo.MatrixSpace(domain, m, n)
+    A = zero(S)
+    for (p_idx, poly) in enumerate(system)
+        for (l_idx, label) in enumerate(labels)
+            A[p_idx, l_idx] = coeff(poly, label)
+        end
+    end
+    kerdim, ker = Nemo.left_kernel(A)
+    relations = [ker[i, :] for i in 1:kerdim]
+    gens(poly_ring), kerdim, ker, relations
+end
diff --git a/src/parametrization.jl b/src/parametrization.jl
index d025c63..cc58e47 100644
--- a/src/parametrization.jl
+++ b/src/parametrization.jl
@@ -55,7 +55,7 @@ function perform_change_of_variables(system, invariants; new_vars_name="y")
 
     newsystem = [Nemo.evaluate(p, substitutions) for p in newsystem]
 
-    return (transform, newsystem)
+    return (Matrix(transform_matrix), transform, newsystem)
 end
 
 """