Added the ability to select a method at runtime if so desired

The Tolerances trait: the a_tol, a_tol_full and r_tol setters move off the
inherent impls onto a trait, with absolute_tolerance and relative_tolerance
getters, so code generic over the integrator can set and read them.

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.4.0"
+version = "0.4.1"
 edition = "2021"
 authors = ["Connor Johnstone"]
 description = "A library for solving differential equations based on the DifferentialEquations.jl julia library."
@@ -24,6 +24,7 @@ roots = "0.0.8"
 approx = "0.5"
 num-dual = "0.15"
 criterion = "0.7.0"
+serde_json = "1.0"
 
 [[bench]]
 name = "simple_1d"

readme.md

@@ -29,6 +29,7 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy
 - **Lazy computation** of extra interpolation stages (Vern7)
 - **PID step size control** as an alternative to the PI controller
 - **Automatic Differentiation**, the same problem can be solved with floats or also differentiated with dual numbers
+- **Runtime integrator selection**, the method can be a value, named in a configuration file rather than fixed as a type
 
 ### Future Improvements
 
@@ -83,6 +84,36 @@ let solution = problem.solve();
 let _interpolated_answer = solution.interpolate(4.4);
 ```
 
+## Choosing the integrator at runtime
+
+The integrator is normally a type parameter, fixed when the code is compiled. `Method` is an enum
+over the integrators here that itself implements `Integrator`, so the choice can be a value instead:
+
+```rust
+use nalgebra::Vector1;
+use ordinary_diffeq::prelude::*;
+
+fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+    -y
+}
+
+// `MethodName` derives `Serialize`/`Deserialize`, so this can come out of a
+// configuration file instead of being written here.
+let name: MethodName = serde_json::from_str(r#""Vern7""#).unwrap();
+let method: Method<1> = name.instantiate().a_tol(1e-12).r_tol(1e-12);
+
+let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(1.0), ());
+let mut problem = Problem::new(ode, method, PIController::default());
+let solution = problem.solve();
+
+// And the method that ran can be recorded in a form that selects it again
+assert_eq!(method.name(), MethodName::Vern7);
+```
+
+Dispatch is a `match` on the variant, not a trait object, so nothing is erased: the same `Method`
+works at a dual scalar and differentiates exactly as the integrator it wraps does. A `Vec<Method<D>>`
+holds several methods for comparison over the same problem.
+
 ## Automatic differentiation
 
 State, time and step size are generic over the `Scalar` trait. `f64` implements it by default,

roadmap/README.md

@@ -68,12 +68,14 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Effort**: Small
   - **Status**: All four integrators implement `Tolerances`; no call site changed
 
-- [ ] **[Runtime Integrator Selection](features/41-runtime-integrator-selection.md)**
+- [x] **[Runtime Integrator Selection](features/41-runtime-integrator-selection.md)** ✅ COMPLETED
   - Structural facts move from associated constants to `&self` methods, so an
     enum can dispatch and an integrator can be chosen at runtime or named in a
     config file
   - **Dependencies**: None (pairs with Tolerances Trait)
   - **Effort**: Medium
+  - **Status**: `Method` and `MethodName` added; benches measured at -0.7% to
+    +0.4%, within criterion's noise threshold
 
 ### Callbacks
 
@@ -354,11 +356,11 @@ Each algorithm implementation should include:
 ## Progress Tracking
 
 Total Features: 41
-- Tier 1: 10 features (5/10 complete) ✅
+- Tier 1: 10 features (6/10 complete) ✅
 - Tier 2: 11 features (0/11 complete)
 - Tier 3: 20 features (0/20 complete)
 
-**Overall Progress: 12.2% (5/41 features complete)**
+**Overall Progress: 14.6% (6/41 features complete)**
 
 ### Completed Features
 1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23)
@@ -366,5 +368,6 @@ Total Features: 41
 3. ✅ PID Controller - Tier 1 (2025-10-24)
 4. ✅ Rosenbrock23 - Tier 1 (2025-10-24)
 5. ✅ Tolerances Trait - Tier 1 (2026-09-30)
+6. ✅ Runtime Integrator Selection - Tier 1 (2026-09-30)
 
 Last updated: 2026-09-30

roadmap/features/41-runtime-integrator-selection.md

@@ -1,6 +1,10 @@
 # Feature: Runtime Integrator Selection
 
-**Status**: proposed (2026-09-30)
+**Status**: completed (2026-09-30)
+
+Measured on the existing criterion benches, constants to methods: `bench_orbit`
+-0.73%, `bench_simple_1d` +0.42%, `bench_interpolation_1d` -0.09%. Criterion
+calls all three within its noise threshold, so the inlining claim below holds.
 
 ## Overview
 

src/integrator/bs3.rs

@@ -4,6 +4,9 @@ use super::super::ode::ODE;
 use super::{Integrator, StepResult, Tolerances};
 use crate::scalar::Scalar;
 
+/// The number of stages BS3 evaluates per step.
+const STAGES: usize = 4;
+
 /// Bogacki-Shampine 3/2 integrator trait for tableau coefficients
 pub trait BS3Integrator<'a> {
     const A: &'a [f64];
@@ -155,14 +158,25 @@ impl<'a, const D: usize> Integrator<D> for BS3<D>
 where
     BS3<D>: BS3Integrator<'a>,
 {
-    const ORDER: usize = 3;
-    const STAGES: usize = 4;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        3
+    }
+
+    fn stages(&self) -> usize {
+        STAGES
+    }
+
+    fn adaptive(&self) -> bool {
+        true
+    }
+
+    fn dense(&self) -> bool {
+        true
+    }
 
     fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
         // Allocate storage for the 4 stages
-        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); STAGES];
 
         // Stage 1: k1 = f(t, y)
         k[0] = (ode.f)(ode.t, ode.y, &ode.params);
@@ -271,13 +285,12 @@ mod tests {
     use nalgebra::Vector1;
 
     #[test]
-    #[allow(clippy::assertions_on_constants)]
     fn test_bs3_creation() {
-        let _bs3: BS3<1> = BS3::new();
-        assert_eq!(BS3::<1>::ORDER, 3);
-        assert_eq!(BS3::<1>::STAGES, 4);
-        assert!(BS3::<1>::ADAPTIVE);
-        assert!(BS3::<1>::DENSE);
+        let bs3: BS3<1> = BS3::new();
+        assert_eq!(bs3.order(), 3);
+        assert_eq!(bs3.stages(), 4);
+        assert!(bs3.adaptive());
+        assert!(bs3.dense());
     }
 
     #[test]

src/integrator/dormand_prince.rs

@@ -4,6 +4,10 @@ use super::{Integrator, StepResult, Tolerances};
 use crate::ode::ODE;
 use crate::scalar::Scalar;
 
+/// The number of stages Dormand-Prince 5(4) evaluates per step. Also the offset
+/// between the 5th and 4th order weights in `B`.
+const STAGES: usize = 7;
+
 /// Integrator Trait
 pub trait DormandPrinceIntegrator<'a> {
     const A: &'a [f64];
@@ -174,24 +178,35 @@ impl<'a, const D: usize> Integrator<D> for DormandPrince45<D>
 where
     DormandPrince45<D>: DormandPrinceIntegrator<'a>,
 {
-    const ORDER: usize = 5;
-    const STAGES: usize = 7;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        5
+    }
+
+    fn stages(&self) -> usize {
+        STAGES
+    }
+
+    fn adaptive(&self) -> bool {
+        true
+    }
+
+    fn dense(&self) -> bool {
+        true
+    }
 
     fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
-        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); STAGES];
         let mut next_y = ode.y;
         let mut err = SVector::<T, D>::zeros();
         let mut rcont5 = SVector::<T, D>::zeros();
         // Do the first of the summations
         k[0] = (ode.f)(ode.t, ode.y, &ode.params);
         next_y += k[0] * T::lift(Self::B[0]) * h;
-        err += k[0] * T::lift(Self::B[0] - Self::B[Self::STAGES]) * h;
+        err += k[0] * T::lift(Self::B[0] - Self::B[STAGES]) * h;
         let rcont1 = ode.y;
         rcont5 += k[0] * h * T::lift(Self::D[0]);
         // Then the rest
-        for i in 1..Self::STAGES {
+        for i in 1..STAGES {
             // Compute the ks
             let mut y_term = SVector::<T, D>::zeros();
             for (j, item) in k.iter().enumerate().take(i) {
@@ -201,12 +216,12 @@ where
 
             // Use that and bis to calculate the y and error terms
             next_y += k[i] * h * T::lift(Self::B[i]);
-            err += k[i] * T::lift(Self::B[i] - Self::B[i + Self::STAGES]) * h;
+            err += k[i] * T::lift(Self::B[i] - Self::B[i + STAGES]) * h;
             rcont5 += k[i] * h * T::lift(Self::D[i]);
         }
         let rcont2 = next_y - ode.y;
         let rcont3 = k[0] * h - rcont2;
-        let rcont4 = rcont2 - k[Self::STAGES - 1] * h - rcont3;
+        let rcont4 = rcont2 - k[STAGES - 1] * h - rcont3;
         let rcont = vec![rcont1, rcont2, rcont3, rcont4, rcont5];
         // Just take the real component of the error
         let mut error_norm_sq = 0.0;

src/integrator/method.rs

@@ -0,0 +1,444 @@
+//! An integrator chosen at runtime.
+//!
+//! Which integrator solves a problem is normally a type parameter, fixed when
+//! the code is compiled. [`Method`] is an enum over the integrators this crate
+//! provides that itself implements [`Integrator`], so a value can stand where a
+//! type would otherwise have to: an accuracy preset can carry a method as well
+//! as a tolerance, a method can be named in a configuration file through
+//! [`MethodName`], and a `Vec<Method<D>>` can hold several for comparison.
+//!
+//! Dispatch is a `match` on the variant, which the compiler resolves per call
+//! site. There is no trait object and no scalar erasure, so autodiff keeps
+//! working through a `Method` exactly as through the integrator it wraps.
+
+use nalgebra::SVector;
+use serde::{Deserialize, Serialize};
+
+use super::bs3::BS3;
+use super::dormand_prince::DormandPrince45;
+use super::rosenbrock::Rosenbrock23;
+use super::vern7::Vern7;
+use super::{Integrator, StepResult, Tolerances};
+use crate::ode::ODE;
+use crate::scalar::Scalar;
+
+/// An integrator chosen at runtime.
+///
+/// Every variant wraps a configured integrator, so tolerances are set on a
+/// `Method` the same way they are set on the integrator itself.
+///
+/// # Example
+///
+/// ```rust
+/// use ordinary_diffeq::prelude::*;
+/// use nalgebra::Vector1;
+///
+/// fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+///     -y
+/// }
+///
+/// // The method is a value here, so it can come from a configuration file, a
+/// // command line flag, or an accuracy preset.
+/// let method: Method<1> = MethodName::Vern7.instantiate().a_tol(1e-12).r_tol(1e-12);
+///
+/// let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(1.0), ());
+/// let mut problem = Problem::new(ode, method, PIController::default());
+/// let solution = problem.solve();
+///
+/// let y = solution.interpolate(1.3)[0];
+/// assert!((y - (-1.3f64).exp()).abs() < 1e-10);
+/// ```
+#[non_exhaustive]
+#[derive(Debug, Clone, Copy)]
+pub enum Method<const D: usize> {
+    /// Bogacki-Shampine 3(2).
+    BS3(BS3<D>),
+    /// Dormand-Prince 5(4).
+    DormandPrince45(DormandPrince45<D>),
+    /// Verner's "most efficient" 7(6).
+    Vern7(Vern7<D>),
+    /// Rosenbrock-W 2(3), for stiff problems.
+    Rosenbrock23(Rosenbrock23<D>),
+}
+
+/// Forwards a `&self` call to whichever integrator the method holds.
+macro_rules! dispatch {
+    ($method:expr, $inner:ident => $call:expr) => {
+        match $method {
+            Method::BS3($inner) => $call,
+            Method::DormandPrince45($inner) => $call,
+            Method::Vern7($inner) => $call,
+            Method::Rosenbrock23($inner) => $call,
+        }
+    };
+}
+
+/// Forwards a consuming call and puts the result back in the same variant.
+macro_rules! dispatch_rewrap {
+    ($method:expr, $inner:ident => $call:expr) => {
+        match $method {
+            Method::BS3($inner) => Method::BS3($call),
+            Method::DormandPrince45($inner) => Method::DormandPrince45($call),
+            Method::Vern7($inner) => Method::Vern7($call),
+            Method::Rosenbrock23($inner) => Method::Rosenbrock23($call),
+        }
+    };
+}
+
+impl<const D: usize> Method<D> {
+    /// Which method this is, without its configuration.
+    ///
+    /// The counterpart to [`MethodName::instantiate`], so a run can record the
+    /// method it used in a form that selects the same one again.
+    pub fn name(&self) -> MethodName {
+        match self {
+            Method::BS3(_) => MethodName::BS3,
+            Method::DormandPrince45(_) => MethodName::DormandPrince45,
+            Method::Vern7(_) => MethodName::Vern7,
+            Method::Rosenbrock23(_) => MethodName::Rosenbrock23,
+        }
+    }
+}
+
+impl<const D: usize> Integrator<D> for Method<D> {
+    fn order(&self) -> usize {
+        dispatch!(self, m => m.order())
+    }
+
+    fn stages(&self) -> usize {
+        dispatch!(self, m => m.stages())
+    }
+
+    fn adaptive(&self) -> bool {
+        dispatch!(self, m => m.adaptive())
+    }
+
+    fn dense(&self) -> bool {
+        dispatch!(self, m => m.dense())
+    }
+
+    fn main_stages(&self) -> usize {
+        dispatch!(self, m => m.main_stages())
+    }
+
+    fn extra_stages(&self) -> usize {
+        dispatch!(self, m => m.extra_stages())
+    }
+
+    fn total_dense_stages(&self) -> usize {
+        dispatch!(self, m => m.total_dense_stages())
+    }
+
+    fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
+        dispatch!(self, m => m.step(ode, h))
+    }
+
+    fn interpolate<T: Scalar>(
+        &self,
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D> {
+        dispatch!(self, m => m.interpolate(t_start, t_end, dense, t))
+    }
+
+    fn compute_extra_stages<T: Scalar, P>(
+        &self,
+        ode: &ODE<T, D, P>,
+        t_start: T,
+        y_start: SVector<T, D>,
+        h: T,
+        main_stages: &[SVector<T, D>],
+    ) -> Vec<SVector<T, D>> {
+        dispatch!(self, m => m.compute_extra_stages(ode, t_start, y_start, h, main_stages))
+    }
+}
+
+impl<const D: usize> Tolerances<D> for Method<D> {
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        dispatch!(self, m => m.absolute_tolerance())
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        dispatch!(self, m => m.relative_tolerance())
+    }
+
+    fn a_tol(self, a_tol: f64) -> Self {
+        dispatch_rewrap!(self, m => m.a_tol(a_tol))
+    }
+
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
+        dispatch_rewrap!(self, m => m.a_tol_full(a_tol))
+    }
+
+    fn r_tol(self, r_tol: f64) -> Self {
+        dispatch_rewrap!(self, m => m.r_tol(r_tol))
+    }
+}
+
+impl<const D: usize> From<BS3<D>> for Method<D> {
+    fn from(integrator: BS3<D>) -> Self {
+        Method::BS3(integrator)
+    }
+}
+
+impl<const D: usize> From<DormandPrince45<D>> for Method<D> {
+    fn from(integrator: DormandPrince45<D>) -> Self {
+        Method::DormandPrince45(integrator)
+    }
+}
+
+impl<const D: usize> From<Vern7<D>> for Method<D> {
+    fn from(integrator: Vern7<D>) -> Self {
+        Method::Vern7(integrator)
+    }
+}
+
+impl<const D: usize> From<Rosenbrock23<D>> for Method<D> {
+    fn from(integrator: Rosenbrock23<D>) -> Self {
+        Method::Rosenbrock23(integrator)
+    }
+}
+
+/// The name of an integrator, without its configuration.
+///
+/// Serializes as the variant name, so `integrator = "Vern7"` in a configuration
+/// file deserializes to [`MethodName::Vern7`] and
+/// [`instantiate`](MethodName::instantiate) turns it into a [`Method`] at
+/// whatever dimension the problem needs.
+///
+/// # Example
+///
+/// ```rust
+/// use ordinary_diffeq::prelude::*;
+///
+/// let name: MethodName = serde_json::from_str(r#""Vern7""#).unwrap();
+/// let method: Method<6> = name.instantiate();
+/// assert_eq!(method.name(), MethodName::Vern7);
+/// ```
+#[non_exhaustive]
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
+pub enum MethodName {
+    /// Bogacki-Shampine 3(2).
+    BS3,
+    /// Dormand-Prince 5(4).
+    DormandPrince45,
+    /// Verner's "most efficient" 7(6).
+    Vern7,
+    /// Rosenbrock-W 2(3), for stiff problems.
+    Rosenbrock23,
+}
+
+impl MethodName {
+    /// Builds the named integrator, at its default tolerances.
+    pub fn instantiate<const D: usize>(&self) -> Method<D> {
+        match self {
+            MethodName::BS3 => Method::BS3(BS3::new()),
+            MethodName::DormandPrince45 => Method::DormandPrince45(DormandPrince45::new()),
+            MethodName::Vern7 => Method::Vern7(Vern7::new()),
+            MethodName::Rosenbrock23 => Method::Rosenbrock23(Rosenbrock23::new()),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::controller::PIController;
+    use crate::problem::Problem;
+    use nalgebra::Vector2;
+
+    /// The harmonic oscillator, which every method solves and which has enough
+    /// curvature that steps land in different places at different tolerances.
+    fn derivative(_t: f64, y: Vector2<f64>, _p: &()) -> Vector2<f64> {
+        Vector2::new(y[1], -y[0])
+    }
+
+    /// Solves to `t = 4` and reports the step times, the states, and the value
+    /// interpolated at ten points that fall strictly inside steps.
+    fn solve<S: Integrator<2> + Copy>(integrator: S) -> (Vec<f64>, Vec<Vector2<f64>>, Vec<f64>) {
+        let ode = crate::ode::ODE::new(&derivative, 0.0, 4.0, Vector2::new(1.0, 0.0), ());
+        let mut problem = Problem::new(ode, integrator, PIController::default());
+        let solution = problem.solve();
+
+        let interpolated = (1..=10)
+            .map(|i| solution.interpolate(i as f64 * 0.37)[0])
+            .collect();
+        (
+            solution.times.clone(),
+            solution.states.clone(),
+            interpolated,
+        )
+    }
+
+    /// Wrapping an integrator in `Method` must not change a single bit of the
+    /// answer, interpolation included — which for Vern7 exercises the lazy
+    /// extra stages, the path `extra_stages` governs.
+    #[test]
+    fn test_method_matches_wrapped_integrator() {
+        macro_rules! assert_identical {
+            ($integrator:expr, $variant:ident) => {{
+                let integrator = $integrator;
+                let (direct_t, direct_y, direct_interp) = solve(integrator);
+                let (wrapped_t, wrapped_y, wrapped_interp) = solve(Method::$variant(integrator));
+
+                assert_eq!(
+                    direct_t,
+                    wrapped_t,
+                    "{} step times differ",
+                    stringify!($variant)
+                );
+                assert_eq!(
+                    direct_y,
+                    wrapped_y,
+                    "{} states differ",
+                    stringify!($variant)
+                );
+                assert_eq!(
+                    direct_interp,
+                    wrapped_interp,
+                    "{} interpolated values differ",
+                    stringify!($variant)
+                );
+                // And the step count is nontrivial, so the comparison means something.
+                assert!(
+                    direct_t.len() > 5,
+                    "{} took too few steps",
+                    stringify!($variant)
+                );
+            }};
+        }
+
+        assert_identical!(BS3::<2>::new().a_tol(1e-8).r_tol(1e-8), BS3);
+        assert_identical!(
+            DormandPrince45::<2>::new().a_tol(1e-10).r_tol(1e-10),
+            DormandPrince45
+        );
+        assert_identical!(Vern7::<2>::new().a_tol(1e-12).r_tol(1e-12), Vern7);
+        assert_identical!(
+            Rosenbrock23::<2>::new().a_tol(1e-6).r_tol(1e-6),
+            Rosenbrock23
+        );
+    }
+
+    /// The structural facts a `Method` reports are the wrapped integrator's.
+    #[test]
+    fn test_method_reports_wrapped_structure() {
+        let vern7 = Vern7::<2>::new();
+        let wrapped = Method::Vern7(vern7);
+        assert_eq!(wrapped.order(), vern7.order());
+        assert_eq!(wrapped.stages(), vern7.stages());
+        assert_eq!(wrapped.main_stages(), vern7.main_stages());
+        assert_eq!(wrapped.extra_stages(), vern7.extra_stages());
+        assert_eq!(wrapped.total_dense_stages(), vern7.total_dense_stages());
+        assert!(wrapped.adaptive() && wrapped.dense());
+
+        // Dormand-Prince has no extra stages, and the default applies through
+        // the enum as it does directly.
+        let dp45 = Method::DormandPrince45(DormandPrince45::<2>::new());
+        assert_eq!(dp45.extra_stages(), 0);
+        assert_eq!(dp45.total_dense_stages(), dp45.main_stages());
+    }
+
+    /// A name round-trips through serde and instantiates the variant it names.
+    #[test]
+    fn test_method_name_round_trip() {
+        for (name, json) in [
+            (MethodName::BS3, r#""BS3""#),
+            (MethodName::DormandPrince45, r#""DormandPrince45""#),
+            (MethodName::Vern7, r#""Vern7""#),
+            (MethodName::Rosenbrock23, r#""Rosenbrock23""#),
+        ] {
+            assert_eq!(serde_json::to_string(&name).unwrap(), json);
+
+            let parsed: MethodName = serde_json::from_str(json).unwrap();
+            assert_eq!(parsed, name);
+
+            let method: Method<3> = parsed.instantiate();
+            assert_eq!(method.name(), name);
+        }
+    }
+
+    /// Tolerance setters forward through the enum and keep the variant.
+    #[test]
+    fn test_tolerances_through_the_enum() {
+        let method: Method<2> = MethodName::BS3.instantiate().a_tol(1e-11).r_tol(1e-9);
+        assert_eq!(method.name(), MethodName::BS3);
+        assert_eq!(
+            method.absolute_tolerance(),
+            SVector::<f64, 2>::repeat(1e-11)
+        );
+        assert_eq!(method.relative_tolerance(), 1e-9);
+
+        let per = SVector::<f64, 2>::new(1e-6, 1e-12);
+        let method = method.a_tol_full(per);
+        assert_eq!(method.absolute_tolerance(), per);
+    }
+
+    /// The heterogeneous collection the type parameter used to rule out.
+    #[test]
+    fn test_collection_of_methods() {
+        let methods: Vec<Method<2>> = [
+            MethodName::BS3,
+            MethodName::DormandPrince45,
+            MethodName::Vern7,
+            MethodName::Rosenbrock23,
+        ]
+        .iter()
+        .map(|name| name.instantiate().a_tol(1e-8).r_tol(1e-8))
+        .collect();
+
+        for method in methods {
+            let (times, states, _) = solve(method);
+            // cos(4) at the end of the interval, to whatever each method manages.
+            let final_y = states.last().unwrap()[0];
+            assert!(
+                (final_y - 4.0f64.cos()).abs() < 1e-3,
+                "{:?} ended at {final_y} after {} steps",
+                method.name(),
+                times.len()
+            );
+        }
+    }
+
+    /// The module header claims a `Method` differentiates exactly as the
+    /// integrator it wraps does, which is only true if nothing erases the
+    /// scalar. Vern7 is the demanding case: interpolating inside a step calls
+    /// the right-hand side again, in dual arithmetic, after the solve.
+    #[cfg(feature = "autodiff")]
+    #[test]
+    fn test_method_preserves_autodiff() {
+        use num_dual::Dual64;
+
+        fn decay<T: crate::scalar::Scalar>(_t: T, y: nalgebra::Vector1<T>) -> nalgebra::Vector1<T> {
+            -y
+        }
+
+        let y0 = Dual64::from_re(1.0).derivative();
+        let ode = crate::ode::ODE::new(
+            &|t, y, _p: &()| decay(t, y),
+            Dual64::lift(0.0),
+            Dual64::lift(2.0),
+            nalgebra::Vector1::new(y0),
+            (),
+        );
+        let method: Method<1> = MethodName::Vern7.instantiate().a_tol(1e-12).r_tol(1e-12);
+        let mut problem = Problem::new(ode, method, PIController::default());
+        let solution = problem.solve();
+
+        let y = solution.interpolate(Dual64::lift(1.3))[0];
+
+        // y(t) = y0 e^-t, so both the value and dy/dy0 are e^-1.3.
+        assert!((y.re - (-1.3f64).exp()).abs() < 1e-10);
+        assert!((y.eps - (-1.3f64).exp()).abs() < 1e-10);
+    }
+
+    /// `From` gets an already-configured integrator into the enum.
+    #[test]
+    fn test_from_integrator() {
+        let method: Method<2> = Vern7::<2>::new().a_tol(1e-9).into();
+        assert_eq!(method.name(), MethodName::Vern7);
+        assert_eq!(method.absolute_tolerance(), SVector::<f64, 2>::repeat(1e-9));
+    }
+}

src/integrator/mod.rs

@@ -5,6 +5,7 @@ use crate::scalar::Scalar;
 
 pub mod bs3;
 pub mod dormand_prince;
+pub mod method;
 pub mod rosenbrock;
 pub mod vern7;
 
@@ -23,27 +24,48 @@ pub struct StepResult<T: Scalar, const D: usize> {
     /// [`Integrator::interpolate`], or `None` from a method without dense
     /// output.
     ///
-    /// For a method with lazy dense output (`EXTRA_STAGES > 0`) the layout is
-    /// `[y_start, y_end, k_1, ..., k_MAIN_STAGES]`, and the extra stages are
+    /// For a method with lazy dense output (`extra_stages() > 0`) the layout is
+    /// `[y_start, y_end, k_1, ..., k_main_stages]`, and the extra stages are
     /// appended by [`Integrator::compute_extra_stages`] on first use.
     pub dense: Option<Vec<SVector<T, D>>>,
 }
 
 /// Integrator Trait
+///
+/// The structural facts about a method are `&self` methods rather than
+/// associated constants, so that a value can answer them. That is what lets
+/// [`Method`](method::Method) wrap the integrators in an enum and dispatch to
+/// whichever one it holds, and so lets an integrator be chosen at runtime. None
+/// of these appears in type position, and for a concrete integrator each call
+/// inlines to the constant it returns.
 pub trait Integrator<const D: usize> {
-    const ORDER: usize;
-    const STAGES: usize;
-    const ADAPTIVE: bool;
-    const DENSE: bool;
+    /// The order of the method's solution.
+    fn order(&self) -> usize;
 
-    /// Number of main stages stored in dense output (default: same as STAGES)
-    const MAIN_STAGES: usize = Self::STAGES;
+    /// The number of stages the method evaluates per step.
+    fn stages(&self) -> usize;
+
+    /// Whether the method estimates its local error, and so can be stepped
+    /// adaptively.
+    fn adaptive(&self) -> bool;
+
+    /// Whether the method produces dense output.
+    fn dense(&self) -> bool;
+
+    /// Number of main stages stored in dense output (default: same as `stages`)
+    fn main_stages(&self) -> usize {
+        self.stages()
+    }
 
     /// Number of extra stages for full-order dense output (default: 0, no extra stages)
-    const EXTRA_STAGES: usize = 0;
+    fn extra_stages(&self) -> usize {
+        0
+    }
 
     /// Total stages when full dense output is computed
-    const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES;
+    fn total_dense_stages(&self) -> usize {
+        self.main_stages() + self.extra_stages()
+    }
 
     /// Advances `ode` by `h` from its current state, without mutating it.
     fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D>;

src/integrator/rosenbrock.rs

@@ -229,10 +229,21 @@ impl<const D: usize> Default for Rosenbrock23<D> {
 }
 
 impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
-    const ORDER: usize = 2;
-    const STAGES: usize = 2;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        2
+    }
+
+    fn stages(&self) -> usize {
+        2
+    }
+
+    fn adaptive(&self) -> bool {
+        true
+    }
+
+    fn dense(&self) -> bool {
+        true
+    }
 
     fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
         let t = ode.t;

src/integrator/vern7.rs

@@ -4,6 +4,13 @@ use super::super::ode::ODE;
 use super::{Integrator, StepResult, Tolerances};
 use crate::scalar::Scalar;
 
+/// The number of main stages Vern7 evaluates per step.
+const STAGES: usize = 10;
+
+/// The number of further stages Vern7 computes, lazily, for full-order dense
+/// output.
+const EXTRA_STAGES: usize = 6;
+
 /// Verner 7 integrator trait for tableau coefficients
 pub trait Vern7Integrator<'a> {
     const A: &'a [f64]; // Lower triangular A matrix (flattened)
@@ -289,24 +296,40 @@ impl<'a, const D: usize> Integrator<D> for Vern7<D>
 where
     Vern7<D>: Vern7Integrator<'a> + Vern7ExtraStages<'a>,
 {
-    const ORDER: usize = 7;
-    const STAGES: usize = 10;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        7
+    }
+
+    fn stages(&self) -> usize {
+        STAGES
+    }
+
+    fn adaptive(&self) -> bool {
+        true
+    }
+
+    fn dense(&self) -> bool {
+        true
+    }
 
     // Lazy dense output configuration
-    const MAIN_STAGES: usize = 10;
-    const EXTRA_STAGES: usize = 6;
+    fn main_stages(&self) -> usize {
+        STAGES
+    }
+
+    fn extra_stages(&self) -> usize {
+        EXTRA_STAGES
+    }
 
     fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
         // Allocate storage for the 10 stages
-        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); STAGES];
 
         // Stage 1: k[0] = f(t, y)
         k[0] = (ode.f)(ode.t, ode.y, &ode.params);
 
         // Compute remaining stages using the A matrix
-        for i in 1..Self::STAGES {
+        for i in 1..STAGES {
             let mut y_temp = ode.y;
             // A matrix is stored in lower triangular form, row by row
             // Row i has i elements (0-indexed), starting at position i*(i-1)/2
@@ -469,7 +492,7 @@ where
                   k9 * b9_theta) * h;
 
         // If extra stages are available, add their contribution for full 7th order accuracy
-        if dense.len() >= 2 + Self::TOTAL_DENSE_STAGES {
+        if dense.len() >= 2 + STAGES + EXTRA_STAGES {
             // Extra stages are at indices 12-17
             let k11 = &dense[12];
             let k12 = &dense[13];
@@ -563,7 +586,7 @@ where
         let k8 = &main_stages[7];
         let k9 = &main_stages[8];
 
-        let mut extra_stages = Vec::with_capacity(Self::EXTRA_STAGES);
+        let mut extra_stages = Vec::with_capacity(EXTRA_STAGES);
 
         // Stage 11: uses k1, k4-k9 (7 coefficients)
         let mut y11 = y_start;

src/lib.rs

@@ -16,6 +16,7 @@ pub mod prelude {
     pub use super::controller::{PIController, PIDController};
     pub use super::integrator::bs3::BS3;
     pub use super::integrator::dormand_prince::DormandPrince45;
+    pub use super::integrator::method::{Method, MethodName};
     pub use super::integrator::rosenbrock::Rosenbrock23;
     pub use super::integrator::vern7::Vern7;
     pub use super::ode::ODE;

src/problem.rs

@@ -52,7 +52,7 @@ where
                 self.controller
                     .set_next_step_guess(TryStep::NotYetAccepted(self.ode.t_end - self.ode.t));
             }
-            let (mut new_y, mut curr_step, mut dense_option) = if S::ADAPTIVE {
+            let (mut new_y, mut curr_step, mut dense_option) = if self.integrator.adaptive() {
                 // First, we try stepping with the "next step guess" to get the error
                 let mut result = self
                     .integrator
@@ -175,13 +175,13 @@ where
                 // Check if we need to compute extra stages for lazy dense output
                 let dense_cell = &self.dense[end_index - 1];
 
-                if S::EXTRA_STAGES > 0 {
+                if self.integrator.extra_stages() > 0 {
                     let needs_extra = {
                         let borrowed = dense_cell.borrow();
                         // Dense array format: [y0, y1, k1, k2, ..., k_main]
-                        // If we have main stages only: 2 + MAIN_STAGES elements
-                        // If we have all stages: 2 + MAIN_STAGES + EXTRA_STAGES elements
-                        borrowed.len() < 2 + S::TOTAL_DENSE_STAGES
+                        // If we have main stages only: 2 + main_stages() elements
+                        // If we have all stages: 2 + total_dense_stages() elements
+                        borrowed.len() < 2 + self.integrator.total_dense_stages()
                     };
 
                     if needs_extra {
@@ -189,7 +189,7 @@ where
                         let mut dense = dense_cell.borrow_mut();
 
                         // Extract main stages (skip y0 and y1 at indices 0 and 1)
-                        let main_stages = &dense[2..2 + S::MAIN_STAGES];
+                        let main_stages = &dense[2..2 + self.integrator.main_stages()];
 
                         // Compute extra stages lazily
                         let extra_stages = self.integrator.compute_extra_stages(