diff --git a/pyproject.toml b/pyproject.toml
index e6dd2cf..583509d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -13,7 +13,7 @@ readme = "README.md"
 requires-python = ">=3.9"
 license = {text = "BSD-3-Clause"}
 dependencies = [
-    "python-libsbml",
+    "python-libsbml>=5.20.4",
     "sympy",
     "pint",
     "lxml>=4.6.4",
diff --git a/sbmlmath/__init__.py b/sbmlmath/__init__.py
index f702380..ea8111e 100644
--- a/sbmlmath/__init__.py
+++ b/sbmlmath/__init__.py
@@ -11,6 +11,7 @@
 from .species_symbol import SpeciesSymbol
 
 __all__ = [
+    "set_math",
     "SBMLMathMLParser",
     "SBMLMathMLPrinter",
     "SpeciesSymbol",
@@ -75,3 +76,33 @@ def sbml_math_to_sympy(
     )
     mathml = libsbml.writeMathMLToString(ast_node)
     return SBMLMathMLParser().parse_str(mathml)
+
+
+def set_math(
+    element: libsbml.SBase,
+    expr: sp.Expr,
+) -> None:
+    """Set the math expression of an SBML object.
+
+    Args:
+        element:
+            The SBML object to set the math expression for.
+        expr:
+            The sympy expression to set as the math expression.
+    """
+    sbml_document = element.getSBMLDocument()
+    mathml = SBMLMathMLPrinter(
+        sbml_level=sbml_document.getLevel(),
+        sbml_version=sbml_document.getVersion(),
+    ).doprint(expr)
+
+    if ast_node := libsbml.readMathMLFromString(mathml):
+        if element.setMath(ast_node) == libsbml.LIBSBML_OPERATION_SUCCESS:
+            return
+        raise ValueError(
+            f"Error setting math expression:\n{expr}\n"
+            f"{mathml}\n{sbml_document.getErrorLog().toString()}"
+        )
+    raise ValueError(
+        f"Unknown error parsing math expression:\n{expr}\n{mathml}"
+    )
diff --git a/tests/test_misc.py b/tests/test_misc.py
index 8edca72..18b9176 100644
--- a/tests/test_misc.py
+++ b/tests/test_misc.py
@@ -4,6 +4,24 @@
 from sbmlmath import *
 
 
+def test_set_math():
+    """Test setting math on SBML objects."""
+    # create some model
+    doc = libsbml.SBMLDocument(3, 1)
+    model = doc.createModel()
+    s = model.createSpecies()
+    s.setId("s")
+    p = model.createParameter()
+    p.setId("p")
+    ia = model.createInitialAssignment()
+    ia.setSymbol("s")
+
+    # actual test
+    expr = sp.sympify("2 * p^2")
+    set_math(ia, expr)
+    assert expr == sbml_math_to_sympy(ia)
+
+
 def test_large_mathml():
     """
     lxml.etree.iterparse simply truncates long input (>2**15 chars);