diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 4aa78ca7f2..a8763c3044 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -96,8 +96,9 @@ include("graphs.jl")
 export Graph, savegraph, complexgraph
 
 # for creating compounds
-include("compound.jl")
+include("chemistry_functionality.jl")
 export @compound
 export components, iscompound, coefficients
+export balance_reaction
 
 end # module
diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl
new file mode 100644
index 0000000000..96b0118f39
--- /dev/null
+++ b/src/chemistry_functionality.jl
@@ -0,0 +1,197 @@
+struct CompoundSpecies end
+struct CompoundComponents end
+struct CompoundCoefficients end
+
+Symbolics.option_to_metadata_type(::Val{:iscompound}) = CompoundSpecies
+Symbolics.option_to_metadata_type(::Val{:components}) = CompoundComponents
+Symbolics.option_to_metadata_type(::Val{:coefficients}) = CompoundCoefficients
+
+macro compound(species_expr, arr_expr...)
+    # Ensure the species name is a valid expression
+    if !(species_expr isa Expr && species_expr.head == :call)
+        error("Invalid species name in @compound macro")
+    end
+
+    # Parse the species name to extract the species name and argument
+    species_name = species_expr.args[1]
+    species_arg = species_expr.args[2]
+
+    # Construct the expressions that define the species
+    species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0),
+                        Expr(:call, species_name, species_arg))
+
+    # Construct the expression to set the iscompound metadata
+    setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
+                                                                       Catalyst.CompoundSpecies,
+                                                                       true))
+
+    # Ensure the expressions are evaluated in the correct scope by escaping them
+    escaped_species_expr = esc(species_expr)
+    escaped_setmetadata_expr = esc(setmetadata_expr)
+
+    # Construct the array from the remaining arguments
+    arr = Expr(:vect, (arr_expr)...)
+    coeffs = []
+    species = []
+
+    for expr in arr_expr
+        if isa(expr, Expr) && expr.head == :call && expr.args[1] == :*
+            push!(coeffs, expr.args[2])
+            push!(species, expr.args[3])
+        else
+            push!(coeffs, 1)
+            push!(species, expr)
+        end
+    end
+
+    coeffs_expr = Expr(:vect, coeffs...)
+    species_expr = Expr(:vect, species...)
+
+    # Construct the expression to set the components metadata
+    setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
+                                                                         Catalyst.CompoundComponents,
+                                                                         $species_expr))
+
+    # Ensure the expression is evaluated in the correct scope by escaping it
+    escaped_setcomponents_expr = esc(setcomponents_expr)
+
+    # Construct the expression to set the coefficients metadata
+    setcoefficients_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
+                                                                           Catalyst.CompoundCoefficients,
+                                                                           $coeffs_expr))
+
+    escaped_setcoefficients_expr = esc(setcoefficients_expr)
+
+    # Return a block that contains the escaped expressions
+    return Expr(:block, escaped_species_expr, escaped_setmetadata_expr,
+                escaped_setcomponents_expr, escaped_setcoefficients_expr)
+end
+
+# Check if a species is a compound
+iscompound(s::Num) = iscompound(MT.value(s))
+function iscompound(s)
+    MT.getmetadata(s, CompoundSpecies, false)
+end
+
+coefficients(s::Num) = coefficients(MT.value(s))
+function coefficients(s)
+    MT.getmetadata(s, CompoundCoefficients)
+end
+
+components(s::Num) = components(MT.value(s))
+function components(s)
+    MT.getmetadata(s, CompoundComponents)
+end
+
+component_coefficients(s::Num) = component_coefficients(MT.value(s))
+function component_coefficients(s)
+    return [c => co for (c, co) in zip(components(s), coefficients(s))]
+end
+
+
+### Balancing Code
+const COMPOUND_OF_COMPOUND_ERROR = ErrorException("Reaction balancing does not currently work for reactions involving compounds of compounds.")
+
+# note this does not correctly handle compounds of compounds currently
+function create_matrix(reaction::Catalyst.Reaction)
+    @unpack substrates, products = reaction
+    unique_atoms = [] # Array to store unique atoms
+    n_atoms = 0
+    ncompounds = length(substrates) + length(products)
+    A = zeros(Int, 0, ncompounds)
+
+    coeffsign = 1
+    jbase = 0
+    for compounds in (substrates, products)
+        for (j2, compound) in enumerate(compounds)
+            j = jbase + j2
+
+            if iscompound(compound)
+                atoms = components(compound)
+                any(iscompound, atoms) && throw(COMPOUND_OF_COMPOUND_ERROR)
+                coeffs = coefficients(compound)
+                (atoms == nothing || coeffs == nothing) && continue
+            else
+                # If not a compound, assume coefficient of 1
+                atoms = [compound]
+                coeffs = [1]
+            end
+
+            for (atom,coeff) in zip(atoms, coeffs)
+                # Extract atom and coefficient from the pair
+                i = findfirst(x -> isequal(x, atom), unique_atoms)
+                if i === nothing
+                    # Add the atom to the atoms array if it's not already present
+                    push!(unique_atoms, atom)
+                    n_atoms += 1
+                    A = [A; zeros(Int, 1, ncompounds)]
+                    i = n_atoms
+                end
+
+                # Adjust coefficient based on whether the compound is a product or substrate
+                coeff *= coeffsign
+
+                A[i, j] = coeff
+            end
+        end
+
+        # update for iterating through products
+        coeffsign = -1
+        jbase = length(substrates)
+    end
+
+    return A
+end
+
+function get_balanced_stoich(reaction::Reaction)
+    # Create the reaction matrix A that is m atoms by n compounds
+    A = create_matrix(reaction)
+
+    # get an integer nullspace basis
+    X = ModelingToolkit.nullspace(A)
+    nullity = size(X, 2)
+
+    stoichvecs = Vector{Vector{Int64}}()
+    for (j, col) in enumerate(eachcol(X))
+        signs = unique(col)
+
+        # If there is only one basis vector and the signs are not all the same this means
+        # we have a solution that would require moving at least one substrate to be a
+        # product (or vice-versa). We therefore do not return anything in this case.
+        # If there are multiple basis vectors we don't currently determine if we can
+        # construct a linear combination giving a solution, so we just return them.
+        if (nullity > 1) || (all(>=(0), signs) || all(<=(0), signs))
+            coefs = abs.(col)
+            common_divisor = reduce(gcd, coefs)
+            coefs .= div.(coefs, common_divisor)
+            push!(stoichvecs, coefs)
+        end
+    end
+
+    return stoichvecs
+end
+
+function balance_reaction(reaction::Reaction)
+    # Calculate the stoichiometric coefficients for the balanced reaction.
+    stoichiometries = get_balanced_stoich(reaction)
+
+    balancedrxs = Vector{Reaction}(undef, length(stoichiometries))
+
+    # Iterate over each stoichiometry vector and create a reaction
+    for (i,stoich) in enumerate(stoichiometries)
+        # Divide the stoichiometry vector into substrate and product stoichiometries.
+        substoich = stoich[1:length(reaction.substrates)]
+        prodstoich = stoich[(length(reaction.substrates) + 1):end]
+
+        # Create a new reaction with the balanced stoichiometries
+        balancedrx = Reaction(reaction.rate, reaction.substrates,
+                              reaction.products, substoich, prodstoich)
+
+        # Add the reaction to the vector of all reactions
+        balancedrxs[i] = balancedrx
+    end
+
+    isempty(balancedrxs) && (@warn "Unable to balance reaction.")
+    (length(balancedrxs) > 1) && (@warn "Infinite balanced reactions from ($reaction) are possible, returning a basis for them. Note that we do not check if they preserve the set of substrates and products from the original reaction.")
+    return balancedrxs
+end
diff --git a/src/compound.jl b/src/compound.jl
deleted file mode 100644
index 01c9081232..0000000000
--- a/src/compound.jl
+++ /dev/null
@@ -1,89 +0,0 @@
-struct CompoundSpecies end
-struct CompoundComponents end
-struct CompoundCoefficients end
-
-Symbolics.option_to_metadata_type(::Val{:iscompound}) = CompoundSpecies
-Symbolics.option_to_metadata_type(::Val{:components}) = CompoundComponents
-Symbolics.option_to_metadata_type(::Val{:coefficients}) = CompoundCoefficients
-
-macro compound(species_expr, arr_expr...)
-    # Ensure the species name is a valid expression
-    if !(species_expr isa Expr && species_expr.head == :call)
-        error("Invalid species name in @compound macro")
-    end
-
-    # Parse the species name to extract the species name and argument
-    species_name = species_expr.args[1]
-    species_arg = species_expr.args[2]
-
-    # Construct the expressions that define the species
-    species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0),
-                        Expr(:call, species_name, species_arg))
-
-    # Construct the expression to set the iscompound metadata
-    setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
-                                                                       Catalyst.CompoundSpecies,
-                                                                       true))
-
-    # Ensure the expressions are evaluated in the correct scope by escaping them
-    escaped_species_expr = esc(species_expr)
-    escaped_setmetadata_expr = esc(setmetadata_expr)
-
-    # Construct the array from the remaining arguments
-    arr = Expr(:vect, (arr_expr)...)
-    coeffs = []
-    species = []
-
-    for expr in arr_expr
-        if isa(expr, Expr) && expr.head == :call && expr.args[1] == :*
-            push!(coeffs, expr.args[2])
-            push!(species, expr.args[3])
-        else
-            push!(coeffs, 1)
-            push!(species, expr)
-        end
-    end
-
-    coeffs_expr = Expr(:vect, coeffs...)
-    species_expr = Expr(:vect, species...)
-
-    # Construct the expression to set the components metadata
-    setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
-                                                                         Catalyst.CompoundComponents,
-                                                                         $species_expr))
-
-    # Ensure the expression is evaluated in the correct scope by escaping it
-    escaped_setcomponents_expr = esc(setcomponents_expr)
-
-    # Construct the expression to set the coefficients metadata
-    setcoefficients_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
-                                                                           Catalyst.CompoundCoefficients,
-                                                                           $coeffs_expr))
-
-    escaped_setcoefficients_expr = esc(setcoefficients_expr)
-
-    # Return a block that contains the escaped expressions
-    return Expr(:block, escaped_species_expr, escaped_setmetadata_expr,
-                escaped_setcomponents_expr, escaped_setcoefficients_expr)
-end
-
-# Check if a species is a compound
-iscompound(s::Num) = iscompound(MT.value(s))
-function iscompound(s)
-    MT.getmetadata(s, CompoundSpecies, false)
-end
-
-coefficients(s::Num) = coefficients(MT.value(s))
-function coefficients(s)
-    MT.getmetadata(s, CompoundCoefficients)
-end
-
-components(s::Num) = components(MT.value(s))
-function components(s)
-    MT.getmetadata(s, CompoundComponents)
-end
-
-component_coefficients(s::Num) = component_coefficients(MT.value(s))
-function component_coefficients(s)
-    return [c => co for (c, co) in zip(components(s), coefficients(s))]
-end
diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl
index fa3059fdcb..be3afbb035 100644
--- a/test/miscellaneous_tests/compound_macro.jl
+++ b/test/miscellaneous_tests/compound_macro.jl
@@ -139,3 +139,4 @@ let
     @test isequal(components(A2),components(B2))
     @test isequal(coefficients(A2), coefficients(B2))
 end
+
diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl
new file mode 100644
index 0000000000..b858cf7aa9
--- /dev/null
+++ b/test/miscellaneous_tests/reaction_balancing.jl
@@ -0,0 +1,406 @@
+using Catalyst, Test
+
+#Check that balancing works.
+let
+    @variables t
+    @parameters k
+    @species H(t) O(t)
+    @compound H2(t) 2H
+    @compound O2(t) 2O
+    @compound H2O(t) 2H 1O
+
+    rx = Reaction(k,[H2,O2],[H2O])
+
+    @test isequal(Catalyst.create_matrix(rx),[2 0 -2; 0 2 -1;])
+    @test isequal(Catalyst.get_balanced_stoich(rx),[[2, 1, 2]])
+
+    balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+let
+    @variables t
+    @parameters k
+    @species C(t) H(t) O(t)
+    @compound O2(t) 2O
+    @compound CO2(t) 1C 2O
+    @compound H2O(t) 2H 1O
+    @compound C6H12O6(t) 6C 12H 6O
+
+    rx = Reaction(k,[CO2,H2O],[C6H12O6,O2])
+
+    @test isequal(Catalyst.create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0;])
+    @test isequal(Catalyst.get_balanced_stoich(rx),[[6, 6, 1, 6]])
+
+    balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1, 6])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, H2O --> H2O
+let
+    @variables t
+    @species H(t) O(t)
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0, [H2O], [H2O], [2], [2])
+
+    balanced_rx = Reaction(1.0, [H2O], [H2O], [1], [1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2H + 1O --> H2O
+let
+    @variables t
+    @species H(t) O(t)
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7])
+
+    balanced_rx = Reaction(1.0, [H,O], [H2O], [2, 1], [1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 1CH4 + 2O2 --> 1CO2 + 2H2O
+let
+    @variables t
+    @species H(t) O(t) C(t)
+    @compound CH4(t) C 4H
+    @compound O2(t) 2O
+    @compound CO2(t) C 2O
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0, [CH4, O2], [CO2, H2O])
+
+    balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, N2 + 3H2 --> 2NH3
+let
+    @variables t
+    @species H(t) N(t)
+    @compound N2(t) 2N
+    @compound H2(t) 2H
+    @compound NH3(t) N 3H
+
+    rx = Reaction(1.0, [N2, H2], [NH3])
+
+    balanced_rx = Reaction(1.0, [N2, H2], [NH3], [1 ,3], [2])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, C2H5OH + CH3COOH --> C4H8O2 + H2O
+let
+    @variables t
+    @species C(t) H(t) O(t)
+    @compound C2H5OH(t) 2C 6H O
+    @compound CH3COOH(t) 2C 4H 2O
+    @compound C4H8O2(t) 4C 8H 2O
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O])
+
+    balanced_rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O], [1, 1], [1, 1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2Ca3PO42 --> 6CaO + 1P4O10
+let
+    @variables t
+    @species Ca(t) P(t) O(t)
+    @compound Ca3PO42(t) 3Ca 2P 8O
+    @compound CaO(t) Ca O
+    @compound P4O10(t) 4P 10O
+
+    rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10])
+
+    balanced_rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10], [2], [6, 1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 4Fe + 3O2 + 6H2O --> 4FeOH3
+let
+    @variables t
+    @species Fe(t) O(t) H(t)
+    @compound O2(t) 2O
+    @compound H2O(t) 2H O
+    @compound FeOH3(t) Fe 3H 3O
+
+    rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3])
+
+    balanced_rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3], [4, 3, 6], [4])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2NaOH + H2SO4 --> Na2SO4 + 2H2O
+let
+    @variables t
+    @species Na(t) O(t) H(t) S(t)
+    @compound SO4(t) S 4O
+    @compound NaOH(t) Na O H
+    @compound H2SO4(t) 2H 1S 4O
+    @compound Na2SO4(t) 2Na 1S 4O
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O])
+
+    balanced_rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O], [2, 1], [1, 2])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2NO2 --> 1N2O4
+let
+    @variables t
+    @species N(t) O(t)
+    @compound NO2(t) N 2O
+    @compound N2O4(t) 2N 4O
+
+    rx = Reaction(1.0, [NO2], [N2O4])
+
+    balanced_rx = Reaction(1.0, [NO2], [N2O4], [2], [1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 1CaCO3 + 2HCl --> CaCl2  + H2O + CO2
+let
+    @variables t
+    @species C(t) H(t) O(t) Ca(t) Cl(t)
+    @compound H2O(t) 2H 1O
+    @compound CO2(t) 1C 2O
+    @compound CaCO3(t) 1Ca 1C 3O
+    @compound HCl(t) 1H 1Cl
+    @compound CaCl2(t) 1Ca 2Cl
+
+    rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O])
+    balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, SiCl4 + 4H2O → H4SiO4 + 4HCl
+let
+    @variables t
+    @species Si(t) Cl(t) H(t) O(t)
+    @compound SiCl4(t) 1Si 4Cl
+    @compound H2O(t) 2H O
+    @compound H4SiO4(t) 4H Si 4O
+    @compound HCl(t) H Cl
+
+    rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl])
+    balanced_rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl], [1,4], [1,4])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2Al + 6HCl → 2AlCl3 + 3H2
+let
+    @variables t
+    @species Al(t) Cl(t) H(t)
+    @compound HCl(t) H Cl
+    @compound AlCl3(t) Al 3Cl
+    @compound H2(t) 2H
+
+    rx = Reaction(1.0,[Al,HCl],[AlCl3,H2])
+    balanced_rx = Reaction(1.0,[Al,HCl],[AlCl3,H2],[2,6], [2,3])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, Na2CO3 + 2HCl → 2NaCl + H2O + CO2
+let
+    @variables t
+    @species Na(t) C(t) O(t) H(t) Cl(t)
+    @compound Na2CO3(t) 2Na C 3O
+    @compound HCl(t) H Cl
+    @compound NaCl(t) Na Cl
+    @compound H2O(t) 2H O
+    @compound CO2(t) C 2O
+
+    rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2])
+    balanced_rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2], [1,2], [2,1,1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2C7H6O2 + 15O2 → 14CO2 + 6H2O
+let
+    @variables t
+    @species C(t) H(t) O(t)
+    @compound C7H6O2(t) 7C 6H 2O
+    @compound O2(t) 2O
+    @compound CO2(t) C 2O
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O])
+    balanced_rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O], [2,15], [14,6])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k,Fe2(SO4)3 + 6KOH → 3K2SO4 + 2Fe(OH)3
+let
+    @variables t
+    @species Fe(t) S(t) O(t) H(t) K(t)
+    @compound Fe2S3O12(t) 2Fe 3S 12O
+    @compound KOH(t) K O H
+    @compound K2SO4(t) 2K S 4O
+    @compound FeO3H3(t) Fe 3O 3H
+
+    rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3]) #5x4 matrix
+    balanced_rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3], [1,6], [3,2])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 2Ca3(PO4)2 + 6SiO2 → P4O10 + 6CaSiO3
+let
+    @variables t
+    @species Ca(t) P(t) O(t) Si(t)
+    @compound Ca3P2O8(t) 3Ca 2P 8O
+    @compound SiO2(t) Si 2O
+    @compound P4O10(t) 4P 10O
+    @compound CaSiO3(t) Ca Si 3O
+
+    rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3]) #5x4 matrix
+    balanced_rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3], [2,6] , [1,6])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, 4KClO3 → 3KClO4 + KCl
+let
+    @variables t
+    @species K(t) Cl(t) O(t)
+    @compound KClO3(t) K Cl 3O
+    @compound KClO4(t) K Cl  4O
+    @compound KCl(t) K Cl
+
+    rx = Reaction(1.0,[KClO3],[KClO4,KCl])
+    balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, Al2(SO4)3 + 3Ca(OH)2 → 2Al(OH)3 + 3CaSO4
+let
+    @variables t
+    @species Al(t) S(t) O(t) Ca(t) O(t) (H)
+    @compound Al2S3O12(t) 2Al 3S 12O
+    @compound CaO2H2(t) Ca 2O 2H
+    @compound AlO3H3(t) Al 3O 3H
+    @compound CaSO4(t) Ca S 4O
+
+    rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4])
+    balanced_rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4], [1,3], [2,3])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, H2SO4 + 8HI → H2S + 4I2 + 4H2O
+let
+    @variables t
+    @species H(t) S(t) O(t) I(t)
+    @compound H2SO4(t) 2H S 4O
+    @compound HI(t) H I
+    @compound H2S(t) 2H S
+    @compound I2(t) 2I
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O])
+    balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# @reaction k, C2H4 + 3O2 ↔ 2CO2 + 2H2O
+let
+    @variables t
+    @species C(t) H(t) O(t)
+    @compound C2H4(t) 2C 4H
+    @compound O2(t) 2O
+    @compound CO2(t) C 2O
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0,[C2H4,O2],[CO2,H2O])
+    balanced_rx = Reaction(1.0,[C2H4,O2],[CO2,H2O],[1,3],[2,2])
+    brxs = balance_reaction(rx)
+    @test length(brxs) == 1
+    @test isequal(balanced_rx, first(brxs))
+end
+
+# Infinite solutions
+let
+    @variables t
+    @species C(t) H(t) O(t)
+    @compound CO(t) C O
+    @compound CO2(t) C 2O
+    @compound H2(t) 2H
+    @compound CH4(t) C 4H
+    @compound H2O(t) 2H O
+
+    rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O])
+    brxs = balance_reaction(rx)
+    @test_logs (:warn, r"Infinite balanced reactions*") match_mode=:any balance_reaction(rx)
+    @test length(brxs) == 2
+end
+
+# No way to balance
+let
+    @variables t
+    @species Fe(t) S(t) O(t) H(t) N(t)
+
+    @compound FeS2(t) Fe 2S
+    @compound HNO3(t) H N 3O
+    @compound Fe2S3O12(t) 2Fe 3S 12O
+    @compound NO(t) N O
+    @compound H2SO4(t) 2H S 4O
+
+    rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4])
+    brxs = balance_reaction(rx)
+    @test_logs (:warn, "Unable to balance reaction.") match_mode=:any balance_reaction(rx)
+    @test isempty(brxs)
+end
+
+# test errors on compounds of compounds
+let
+    @variables t
+    @species C(t) H(t) O(t)
+    @compound CO(t) C O
+    @compound H2(t) 2H
+    @compound COH2(t) CO H2
+
+    rx = Reaction(1.0, [CO, H2], [COH2])
+    @test_throws Catalyst.COMPOUND_OF_COMPOUND_ERROR balance_reaction(rx)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 2bbfa84495..5a772722f2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,6 +23,7 @@ using SafeTestsets
     @time @safetestset "Symbolic Stoichiometry" begin include("miscellaneous_tests/symbolic_stoichiometry.jl") end
     @time @safetestset "Events" begin include("miscellaneous_tests/events.jl") end
     @time @safetestset "Compound species" begin include("miscellaneous_tests/compound_macro.jl") end
+    @time @safetestset "Reaction balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end
     @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end
 
     ### Reaction network analysis. ###