diff --git a/packages/catlog/src/simulate/ode/linear_ode.rs b/packages/catlog/src/simulate/ode/linear_ode.rs
index 01c74785c..12eeed43a 100644
--- a/packages/catlog/src/simulate/ode/linear_ode.rs
+++ b/packages/catlog/src/simulate/ode/linear_ode.rs
@@ -1,44 +1,39 @@
 //! Constant-coefficient linear first-order differential equations.
 
-use nalgebra::{DMatrix, DVector};
+use nalgebra::DMatrix;
+
+use crate::{simulate::ode::NumericalPolynomialSystem, zero::alg::Polynomial};
 
 #[cfg(test)]
 use super::ODEProblem;
-use super::ODESystem;
 
-/// A (constant-coefficient) linear (first-order) dynamical system.
+/// Construct a (constant-coefficient) linear (first-order) dynamical system.
 ///
 /// A system of linear first-order ODEs with constant coefficients; a semantics for
 /// causal loop diagrams.
-#[derive(Clone, Debug, PartialEq)]
-pub struct LinearODESystem {
-    coefficients: DMatrix<f32>,
-}
-
-impl LinearODESystem {
-    /// Constructs a linear ODE system with the given coefficient matrix.
-    pub fn new(A: DMatrix<f32>) -> Self {
-        Self { coefficients: A }
-    }
-}
-
-impl ODESystem for LinearODESystem {
-    fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
-        let A = &self.coefficients;
-        *dx = A * x
+pub fn linear_polynomial_system(coefficients: DMatrix<f32>) -> NumericalPolynomialSystem<u8> {
+    NumericalPolynomialSystem {
+        components: coefficients
+            .row_iter()
+            .map(|row| {
+                row.iter()
+                    .enumerate()
+                    .map(|(j, a)| Polynomial::generator(j) * *a)
+                    .sum::<Polynomial<_, _, _>>()
+            })
+            .collect(),
     }
 }
 
 #[cfg(test)]
-pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem<LinearODESystem> {
+pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem<NumericalPolynomialSystem<u8>> {
     use nalgebra::{dmatrix, dvector};
 
     let A = dmatrix![-0.3,  0.0,  0.0;
                       0.0,  0.0,  0.5;
                       1.0, -2.0,  0.0];
-    let system = LinearODESystem::new(A);
     let initial = dvector![2.0, 1.0, 1.0];
-    ODEProblem::new(system, initial).end_time(10.0)
+    ODEProblem::new(linear_polynomial_system(A), initial).end_time(10.0)
 }
 
 #[cfg(test)]
diff --git a/packages/catlog/src/simulate/ode/lotka_volterra.rs b/packages/catlog/src/simulate/ode/lotka_volterra.rs
index 82561725a..b041ab295 100644
--- a/packages/catlog/src/simulate/ode/lotka_volterra.rs
+++ b/packages/catlog/src/simulate/ode/lotka_volterra.rs
@@ -1,45 +1,45 @@
 //! Lotka-Volterra differential equations.
 
+use crate::simulate::ode::NumericalPolynomialSystem;
+use crate::zero::alg::Polynomial;
 use nalgebra::{DMatrix, DVector};
 
 #[cfg(test)]
 use super::ODEProblem;
-use super::ODESystem;
 
-/// A Lotka-Volterra dynamical system.
+/// Construct a Lotka-Volterra dynamical system.
 ///
 /// A system of ODEs that is affine in its *logarithmic* derivative. These are
 /// sometimes called the "generalized Lotka-Volterra equations." For more, see
 /// [Wikipedia](https://en.wikipedia.org/wiki/Generalized_Lotka%E2%80%93Volterra_equation).
-#[derive(Clone, Debug, PartialEq)]
-pub struct LotkaVolterraSystem {
+pub fn lotka_volterra_system(
     interaction_coeffs: DMatrix<f32>,
     growth_rates: DVector<f32>,
-}
-
-impl LotkaVolterraSystem {
-    /// Constructs a new Lokta-Volterra system with the given parameters.
-    pub fn new(A: DMatrix<f32>, b: DVector<f32>) -> Self {
-        Self { interaction_coeffs: A, growth_rates: b }
-    }
-}
-
-impl ODESystem for LotkaVolterraSystem {
-    fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
-        let A = &self.interaction_coeffs;
-        let b = &self.growth_rates;
-        *dx = (A * x + b).component_mul(x);
+) -> NumericalPolynomialSystem<u8> {
+    NumericalPolynomialSystem {
+        components: interaction_coeffs
+            .row_iter()
+            .enumerate()
+            .zip(growth_rates.into_iter())
+            .map(|((i, row), rate)| {
+                Polynomial::<_, f32, _>::generator(i)
+                    * (row
+                        .iter()
+                        .enumerate()
+                        .map(|(j, a)| Polynomial::generator(j) * *a)
+                        .sum::<Polynomial<_, _, _>>()
+                        + *rate)
+            })
+            .collect(),
     }
 }
 
 #[cfg(test)]
-pub(crate) fn create_predator_prey() -> ODEProblem<LotkaVolterraSystem> {
+pub(crate) fn create_predator_prey() -> ODEProblem<NumericalPolynomialSystem<u8>> {
     let A = DMatrix::from_row_slice(2, 2, &[0.0, -1.0, 1.0, 0.0]);
     let b = DVector::from_column_slice(&[2.0, -1.0]);
-    let sys = LotkaVolterraSystem::new(A, b);
-
     let x0 = DVector::from_column_slice(&[1.0, 1.0]);
-    ODEProblem::new(sys, x0).end_time(10.0)
+    ODEProblem::new(lotka_volterra_system(A, b), x0).end_time(10.0)
 }
 
 #[cfg(test)]
diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs
index 98c5a71dc..518371987 100644
--- a/packages/catlog/src/simulate/ode/polynomial.rs
+++ b/packages/catlog/src/simulate/ode/polynomial.rs
@@ -184,6 +184,7 @@ where
 ///
 /// Such a system is ready for use in numerical solvers: the coefficients are
 /// floating point numbers and the variables are consecutive integer indices.
+#[derive(PartialEq, Debug)]
 pub struct NumericalPolynomialSystem<Exp> {
     /// Components of the vector field.
     pub components: Vec<Polynomial<usize, f32, Exp>>,
diff --git a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs
index 40c06ed17..45798f420 100644
--- a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs
+++ b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs
@@ -14,7 +14,7 @@ use tsify::Tsify;
 
 use super::{ODEAnalysis, SignedCoefficientBuilder};
 use crate::dbl::model::DiscreteDblModel;
-use crate::simulate::ode::{LinearODESystem, ODEProblem};
+use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, linear_polynomial_system};
 use crate::{one::QualifiedPath, zero::QualifiedName};
 
 /// Data defining a linear ODE problem for a model.
@@ -47,7 +47,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
         &self,
         model: &DiscreteDblModel,
         data: LinearODEProblemData,
-    ) -> ODEAnalysis<LinearODESystem> {
+    ) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
         let (matrix, ob_index) = self.build_matrix(model, &data.coefficients);
         let n = ob_index.len();
 
@@ -56,7 +56,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
             .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
         let x0 = DVector::from_iterator(n, initial_values);
 
-        let system = LinearODESystem::new(matrix);
+        let system = linear_polynomial_system(matrix);
         let problem = ODEProblem::new(system, x0).end_time(data.duration);
         ODEAnalysis::new(problem, ob_index)
     }
diff --git a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs
index 323e32701..dd61c5ce2 100644
--- a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs
+++ b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs
@@ -14,7 +14,7 @@ use tsify::Tsify;
 
 use super::{ODEAnalysis, SignedCoefficientBuilder};
 use crate::dbl::model::DiscreteDblModel;
-use crate::simulate::ode::{LotkaVolterraSystem, ODEProblem};
+use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, lotka_volterra_system};
 use crate::{one::QualifiedPath, zero::QualifiedName};
 
 /// Data defining a Lotka-Volterra ODE problem for a model.
@@ -51,7 +51,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
         &self,
         model: &DiscreteDblModel,
         data: LotkaVolterraProblemData,
-    ) -> ODEAnalysis<LotkaVolterraSystem> {
+    ) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
         let (matrix, ob_index) = self.build_matrix(model, &data.interaction_coeffs);
         let n = ob_index.len();
 
@@ -64,8 +64,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
             .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
         let x0 = DVector::from_iterator(n, initial_values);
 
-        let system = LotkaVolterraSystem::new(matrix, b);
-        let problem = ODEProblem::new(system, x0).end_time(data.duration);
+        let problem = ODEProblem::new(lotka_volterra_system(matrix, b), x0).end_time(data.duration);
         ODEAnalysis::new(problem, ob_index)
     }
 }
diff --git a/packages/catlog/src/zero/alg.rs b/packages/catlog/src/zero/alg.rs
index c994653cb..8c422cf72 100644
--- a/packages/catlog/src/zero/alg.rs
+++ b/packages/catlog/src/zero/alg.rs
@@ -1,9 +1,9 @@
 //! Commutative algebra and polynomials.
 
-use num_traits::{One, Pow, Zero};
+use num_traits::{One, Pow, Zero, one, zero};
 use std::collections::BTreeMap;
 use std::fmt::Display;
-use std::iter::{Product, Sum};
+use std::iter::{self, Product, Sum};
 use std::ops::{Add, AddAssign, Mul, Neg};
 
 use derivative::Derivative;
@@ -35,7 +35,7 @@ pub trait CommAlg: CommRing + Module<Ring = Self::R> {
 /// In abstract terms, polynomials with coefficients valued in a [commutative
 /// ring](super::rig::CommRing) *R* are the free [commutative algebra](CommAlg)
 /// over *R*.
-#[derive(Clone, PartialEq, Eq, Derivative)]
+#[derive(Clone, PartialEq, Eq, Derivative, Debug)]
 #[derivative(Default(bound = ""))]
 pub struct Polynomial<Var, Coef, Exp>(Combination<Monomial<Var, Exp>, Coef>);
 
@@ -183,6 +183,29 @@ where
 // XXX: Lots of boilerplate to delegate the module structure of `Polynomial` to
 // `Combination` because these traits cannot be derived automatically.
 
+impl<Var, Coef, Exp> Sum for Polynomial<Var, Coef, Exp>
+where
+    Var: Ord,
+    Coef: AdditiveMonoid,
+    Exp: Ord,
+{
+    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
+        iter.fold(zero(), |acc, x| acc + x)
+    }
+}
+
+impl<Var, Coef, Exp> Add<Coef> for Polynomial<Var, Coef, Exp>
+where
+    Var: Ord,
+    Coef: Add<Output = Coef>,
+    Exp: Ord + Zero,
+{
+    type Output = Polynomial<Var, Coef, Exp>;
+    fn add(self, a: Coef) -> Self::Output {
+        Polynomial(self.0 + iter::once((a, one::<Monomial<_, _>>())).collect())
+    }
+}
+
 impl<Var, Coef, Exp> AddAssign<(Coef, Monomial<Var, Exp>)> for Polynomial<Var, Coef, Exp>
 where
     Var: Ord,
@@ -210,7 +233,7 @@ where
 impl<Var, Coef, Exp> Zero for Polynomial<Var, Coef, Exp>
 where
     Var: Ord,
-    Coef: Add<Output = Coef> + Zero,
+    Coef: Zero,
     Exp: Ord,
 {
     fn zero() -> Self {