Release 0.6.1

Documentation and housekeeping only; the library behaves exactly as 0.6.0 does.
Three rustdoc warnings in the Jacobian comment are gone, so the formulas there
render as code rather than as attempted links, the tree is rustfmt-clean with the
Butcher tableaux and interpolation coefficients marked to keep their hand
alignment, and the git hooks stay out of the package.

Cargo.toml

@@ -1,13 +1,13 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.3.1"
+version = "0.6.1"
 edition = "2021"
 authors = ["Connor Johnstone"]
 description = "A library for solving differential equations based on the DifferentialEquations.jl julia library."
 readme = "readme.md"
 repository = "https://git.rcjohnstone.com/connor/ordinary-diffeq"
 license = "MIT"
-exclude = ["roadmap/"]
+exclude = ["roadmap/", "hooks/"]
 
 [features]
 default = ["autodiff"]
@@ -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"

benches/bs3_vs_dp5.rs

@@ -194,7 +194,11 @@ fn bench_interpolation(c: &mut Criterion) {
             let mut problem = Problem::new(ode, bs3, controller);
             let solution = problem.solve();
             // Interpolate at 100 points
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 
@@ -206,7 +210,11 @@ fn bench_interpolation(c: &mut Criterion) {
             let mut problem = Problem::new(ode, dp45, controller);
             let solution = problem.solve();
             // Interpolate at 100 points
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 

benches/orbit.rs

@@ -1,7 +1,7 @@
 use criterion::{criterion_group, criterion_main, Criterion};
 
-use ordinary_diffeq::prelude::*;
 use nalgebra::Vector6;
+use ordinary_diffeq::prelude::*;
 
 fn bench_orbit(c: &mut Criterion) {
     let mu = 3.98600441500000e14;

benches/simple_1d.rs

@@ -1,7 +1,7 @@
 use criterion::{criterion_group, criterion_main, Criterion};
 
-use ordinary_diffeq::prelude::*;
 use nalgebra::Vector1;
+use ordinary_diffeq::prelude::*;
 
 fn bench_simple_1d(c: &mut Criterion) {
     type Params = (f64,);

benches/vern7_comparison.rs

@@ -172,7 +172,11 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
             let mut problem = Problem::new(ode, vern7, controller);
             let solution = problem.solve();
             // Interpolate at 100 points - first one computes extra stages
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 
@@ -183,7 +187,11 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
             let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
             let mut problem = Problem::new(ode, dp45, controller);
             let solution = problem.solve();
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 

hooks/pre-commit

@@ -0,0 +1,11 @@
+#!/bin/sh
+# Fast checks only. On a warm tree these are about a second between them; a hook
+# that takes a minute is a hook that gets --no-verify. The doctests are in
+# pre-push instead, since each one links its own binary and they cost ~15s.
+#
+# Install: git config core.hooksPath hooks
+set -e
+
+cargo fmt --check
+cargo clippy --all-targets -- -D warnings
+cargo test --quiet --lib

hooks/pre-push

@@ -0,0 +1,9 @@
+#!/bin/sh
+# Everything the pre-commit hook is too slow to run. The doctests matter here
+# because the readme is the crate's front page, so its examples are compiled.
+#
+# Install: git config core.hooksPath hooks
+set -e
+
+cargo test --quiet
+cargo doc --quiet --no-deps

readme.md

@@ -23,12 +23,14 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy
 ### Other Features
 
 - **Adaptive time stepping** with PI controller
-- **Callback events** with zero-crossing detection
+- **Callback events** with zero-crossing detection, whose effects can change the state, change the
+  parameters, or end the solution at the event
 - **Dense output interpolation** at any time point
 - **Parameters** in derivative and callback functions
 - **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
 
@@ -70,19 +72,98 @@ let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
 let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
 let controller = PIController::default();
 
-let value_too_high = FnCallback {
+let mut value_too_high = FnCallback {
     event: &|t: f64, _y: &Vector2<f64>, _p: &Params| { 5.0 - t },
     effect: &stop,
 };
 
 // Solve the problem
-let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high);
+let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut value_too_high);
 let solution = problem.solve();
 
 // Can interpolate solutions to whatever you want
 let _interpolated_answer = solution.interpolate(4.4);
 ```
 
+## Events that change the state
+
+A callback watches a function of the state for a zero crossing. When one is found the step is
+shortened to land on the crossing, the state there is recorded, and the effect runs with the ODE at
+that point — so the effect can move the state, and the integration continues from where it leaves it.
+
+That is what it takes to integrate through a discontinuity. A rotating body parameterized by
+Modified Rodrigues Parameters, for instance, has a singularity at a full turn, avoided by swapping
+the parameter vector for `-sigma / |sigma|^2` — the same attitude, different parameters — whenever
+`|sigma|` passes 1:
+
+```rust
+use nalgebra::{Matrix3, Vector3};
+use ordinary_diffeq::prelude::*;
+
+/// MRP kinematics for a constant angular velocity.
+fn attitude_rate(_t: f64, s: Vector3<f64>, omega: &Vector3<f64>) -> Vector3<f64> {
+    let skew = Matrix3::new(0.0, -s[2], s[1], s[2], 0.0, -s[0], -s[1], s[0], 0.0);
+    let b = Matrix3::identity() * (1.0 - s.dot(&s)) + skew * 2.0 + s * s.transpose() * 2.0;
+    b * omega * 0.25
+}
+
+let mut shadow_set = FnCallback {
+    event: &|_t: f64, s: &Vector3<f64>, _p: &Vector3<f64>| 1.0 - s.norm(),
+    effect: &|ode: &mut ODE<'_, f64, 3, Vector3<f64>>| {
+        ode.y = -ode.y / ode.y.norm_squared();
+    },
+};
+
+let ode = ODE::new(&attitude_rate, 0.0, 5.0, Vector3::zeros(), Vector3::new(0.0, 0.0, 1.0));
+let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+let mut problem = Problem::new(ode, vern7, PIController::default())
+    .with_callback(&mut shadow_set);
+let solution = problem.solve();
+
+// The parameters never approach the singularity, having swapped at half a turn.
+assert!(solution.states.iter().all(|s| s.norm() <= 1.0 + 1e-9));
+```
+
+A jump like that one puts the same time in the solution twice, holding the state on each side of it,
+and makes the step size controller forget a history that described the solution before the jump.
+Interpolating at the jump reports the state the effect left, the one the trajectory goes on with.
+`effect` takes `&mut self`, so a callback that implements the trait itself rather than going through
+`FnCallback` can record that it has fired. `stop` ends the solution at the crossing, with the
+crossing as its last point and the problem's own `t_end` left alone.
+
+If two events fall inside one step the earlier one fires and the step ends there; the later one is
+found again on the step after it.
+
+## 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/GETTING_STARTED.md

@@ -127,6 +127,27 @@ This is a solid foundation! The roadmap builds on this.
 
 ## Development Tips
 
+### Checks, and the Hooks That Run Them
+
+The tree is rustfmt-clean, clippy-clean under `-D warnings`, and builds its docs
+without warnings. All three stay that way, and a pre-existing warning is not an
+exemption — fix it or mark it deliberately, as the Butcher tableaux are marked
+`#[rustfmt::skip]` so rustfmt leaves their hand alignment alone.
+
+The hooks in `hooks/` run them. Git does not track where it looks for hooks, so
+each clone needs one command:
+
+```sh
+git config core.hooksPath hooks
+```
+
+`pre-commit` is `cargo fmt --check`, `cargo clippy --all-targets -- -D warnings`
+and `cargo test --lib`, about a second between them on a warm tree. `pre-push`
+adds the full `cargo test` and `cargo doc --no-deps`, which cost around 25
+seconds — the doctests are slow because each one links its own binary, and they
+matter because the readme is the crate's front page and its examples are
+compiled.
+
 ### Testing Strategy
 
 Every feature should have:

roadmap/README.md

@@ -8,7 +8,8 @@ The library currently has:
 - ✅ Dormand-Prince 4(5) adaptive explicit RK method with dense output
 - ✅ Tsit5 explicit RK method with dense output
 - ✅ PI step size controller
-- ✅ Basic continuous callbacks with zero-crossing detection
+- ✅ Continuous callbacks with zero-crossing detection, whose effects can change
+  the state, the parameters, or end the solution
 - ✅ Solution interpolation interface
 - ✅ Generic over array dimensions at compile-time
 - ✅ Support for ordinary and second-order ODE problems
@@ -58,8 +59,39 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Dependencies**: None
   - **Effort**: Small
 
+### API Ergonomics
+
+- [x] **[Tolerances Trait](features/40-tolerances-trait.md)** ✅ COMPLETED
+  - Move `a_tol` / `a_tol_full` / `r_tol` off the inherent impls into a trait, and
+    add getters, so a generic `S` can set and read them
+  - Raised by `astrokit`, which declared its own copy as a workaround
+  - **Dependencies**: None
+  - **Effort**: Small
+  - **Status**: All four integrators implement `Tolerances`; no call site changed
+
+- [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
 
+- [x] **[Callback Effects That Can Change the State](features/42-callback-effects.md)** ✅ COMPLETED
+  - An effect's write to the state was overwritten by the step it was located in,
+    so a callback could stop a solution but not change one
+  - `effect(&mut self)`, termination at the crossing, both sides of a jump
+    recorded, the controller reset across one, and a defined order for two events
+    in one step
+  - Raised by `astrokit`, which needs a Modified Rodrigues Parameter shadow set
+    swap mid-propagation
+  - **Dependencies**: None
+  - **Effort**: Medium
+  - **Status**: All success criteria met; benches unchanged within run-to-run drift
+
 - [ ] **[Discrete Callbacks](features/05-discrete-callbacks.md)**
   - Event detection based on conditions (not zero-crossings)
   - Useful for time-based events, iteration counts, etc.
@@ -74,13 +106,13 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Solution Interface
 
-- [ ] **[Saveat Functionality](features/07-saveat.md)**
+- [ ] **[Saveat Functionality](features/07-saveat-functionality.md)**
   - Save solution at specific timepoints
   - Dense vs sparse saving strategies
   - **Dependencies**: None
   - **Effort**: Medium
 
-- [ ] **[Solution Derivatives](features/08-solution-derivatives.md)**
+- [ ] **Solution Derivatives** (no plan written yet)
   - Access derivatives at any time point via interpolation
   - `solution.derivative(t)` interface
   - **Dependencies**: None
@@ -92,7 +124,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Algorithms
 
-- [ ] **[DP8 (Dormand-Prince 8th order)](features/09-dp8-method.md)**
+- [ ] **DP8 (Dormand-Prince 8th order)** (no plan written yet)
   - 8th order explicit RK for very high accuracy
   - Complements Vern7 for algorithm selection
   - **Dependencies**: None
@@ -104,7 +136,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Dependencies**: Linear solver infrastructure, Nordsieck representation
   - **Effort**: Large
 
-- [ ] **[Rodas4/Rodas5P](features/11-rodas-methods.md)**
+- [ ] **[Rodas4/Rodas5P](features/11-rodas4-rodas5p-methods.md)**
   - Higher-order Rosenbrock methods (4th/5th order)
   - Better accuracy for stiff problems
   - **Dependencies**: Rosenbrock23
@@ -119,7 +151,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Dependencies**: At least one stiff solver (Rosenbrock23 or FBDF)
   - **Effort**: Large
 
-- [ ] **[Default Algorithm Selection](features/13-default-algorithm.md)**
+- [ ] **[Default Algorithm Selection](features/13-default-algorithm-selection.md)**
   - Smart defaults based on problem characteristics
   - `solve(problem)` without specifying algorithm
   - **Dependencies**: Auto-switching, multiple algorithms
@@ -127,7 +159,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Initialization
 
-- [ ] **[Automatic Initial Step Size](features/14-initial-stepsize.md)**
+- [ ] **[Automatic Initial Step Size](features/14-automatic-initial-step-size.md)**
   - Algorithm to determine good initial dt
   - Based on local Lipschitz estimate
   - **Dependencies**: None
@@ -135,19 +167,19 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Callbacks
 
-- [ ] **[PresetTimeCallback](features/15-preset-time-callback.md)**
+- [ ] **[PresetTimeCallback](features/15-presettimecallback.md)**
   - Trigger callbacks at specific predetermined times
   - Important for time-varying forcing functions
   - **Dependencies**: Discrete callbacks
   - **Effort**: Small
 
-- [ ] **[TerminateSteadyState](features/16-terminate-steady-state.md)**
+- [ ] **[TerminateSteadyState](features/16-terminatesteadystate.md)**
   - Auto-detect when solution reaches steady state
   - Stop integration early
   - **Dependencies**: Discrete callbacks
   - **Effort**: Small
 
-- [ ] **[SavingCallback](features/17-saving-callback.md)**
+- [ ] **[SavingCallback](features/17-savingcallback.md)**
   - Custom saving logic beyond default
   - For memory-efficient large-scale simulations
   - **Dependencies**: CallbackSet
@@ -175,7 +207,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Algorithms
 
-- [ ] **[Low-Storage Runge-Kutta](features/20-low-storage-rk.md)**
+- [ ] **[Low-Storage Runge-Kutta](features/20-low-storage-runge-kutta.md)**
   - 2N, 3N, 4N storage variants
   - Critical for large-scale problems
   - **Dependencies**: None
@@ -193,7 +225,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Dependencies**: Second-order ODE infrastructure (already exists)
   - **Effort**: Medium
 
-- [ ] **[Verner Methods Suite](features/23-verner-methods.md)**
+- [ ] **[Verner Methods Suite](features/23-verner-methods-suite.md)**
   - Vern6, Vern8, Vern9
   - Complete Verner family for different accuracy needs
   - **Dependencies**: Vern7
@@ -239,19 +271,19 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Advanced Callbacks
 
-- [ ] **[VectorContinuousCallback](features/30-vector-continuous-callback.md)**
+- [ ] **[VectorContinuousCallback](features/30-vectorcontinuouscallback.md)**
   - Multiple simultaneous event detection
   - More efficient than multiple callbacks
   - **Dependencies**: CallbackSet
   - **Effort**: Medium
 
-- [ ] **[PositiveDomain](features/31-positive-domain.md)**
+- [ ] **[PositiveDomain](features/31-positivedomain.md)**
   - Enforce positivity constraints
   - Important for physical systems
   - **Dependencies**: CallbackSet
   - **Effort**: Small
 
-- [ ] **[ManifoldProjection](features/32-manifold-projection.md)**
+- [ ] **[ManifoldProjection](features/32-manifoldprojection.md)**
   - Project solution onto constraint manifolds
   - For constrained mechanical systems
   - **Dependencies**: CallbackSet
@@ -259,14 +291,14 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Infrastructure
 
-- [ ] **[Nonlinear Solver Infrastructure](features/33-nonlinear-solver.md)**
+- [ ] **[Nonlinear Solver Infrastructure](features/33-nonlinear-solver-infrastructure.md)**
   - Newton's method
   - Quasi-Newton methods
   - Generic nonlinear solver interface
   - **Dependencies**: Linear solver infrastructure
   - **Effort**: Large
 
-- [ ] **[Krylov Linear Solvers](features/34-krylov-solvers.md)**
+- [ ] **[Krylov Linear Solvers](features/34-krylov-linear-solvers.md)**
   - GMRES, BiCGStab
   - For large sparse systems
   - **Dependencies**: Linear solver infrastructure
@@ -336,17 +368,20 @@ Each algorithm implementation should include:
 
 ## Progress Tracking
 
-Total Features: 38
-- Tier 1: 8 features (4/8 complete) ✅
-- Tier 2: 12 features (0/12 complete)
-- Tier 3: 18 features (0/18 complete)
+Total Features: 42
+- Tier 1: 11 features (7/11 complete) ✅
+- Tier 2: 11 features (0/11 complete)
+- Tier 3: 20 features (0/20 complete)
 
-**Overall Progress: 10.5% (4/38 features complete)**
+**Overall Progress: 16.7% (7/42 features complete)**
 
 ### Completed Features
 1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23)
 2. ✅ Vern7 (Verner 7th order) - Tier 1 (2025-10-24)
 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)
+7. ✅ Callback Effects That Can Change the State - Tier 1 (2026-10-01)
 
-Last updated: 2025-10-24
+Last updated: 2026-10-01

roadmap/features/40-tolerances-trait.md

@@ -0,0 +1,129 @@
+# Feature: Tolerances Trait
+
+**Status**: completed (2026-09-30)
+
+## Overview
+
+Error tolerances are set through inherent builder methods on each integrator
+rather than through a trait:
+
+```rust
+impl<const D: usize> Vern7<D> where Vern7<D>: Integrator<D> {
+    pub fn a_tol(mut self, a_tol: f64) -> Self { .. }
+    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self { .. }
+    pub fn r_tol(mut self, r_tol: f64) -> Self { .. }
+}
+```
+
+`DormandPrince45`, `BS3` and `Rosenbrock23` each have the same four methods,
+written out separately. There is no trait bound that says "this integrator has
+tolerances", so a downstream crate generic over `S: Integrator<D>` cannot set
+them, and cannot read back what they were.
+
+## Why This Feature Matters
+
+A library wrapping this one wants to offer accuracy presets:
+
+```rust
+Propagation::new(..).accuracy(Accuracy::Medium)   // sets a_tol and r_tol
+Propagation::new(..).a_tol(1e-10)
+```
+
+That is impossible over a generic `S`. `astrokit` worked around it by declaring
+its own trait and implementing it for three of the four integrators:
+
+```rust
+pub trait Tolerances<const D: usize>: Integrator<D> {
+    fn with_a_tol(self, a_tol: f64) -> Self;
+    fn with_a_tol_per_component(self, a_tol: SVector<f64, D>) -> Self;
+    fn with_r_tol(self, r_tol: f64) -> Self;
+}
+impl<const D: usize> Tolerances<D> for Vern7<D> { .. }
+```
+
+That works, but the trait belongs here. Downstream it has three problems: any
+integrator added to this crate is invisible to it until the downstream crate is
+updated, a second dependent would declare the same trait again, and the orphan
+rule means only crates that own the trait can extend it.
+
+Reading tolerances back matters separately. A run record that says which
+integrator ran but not at what tolerance does not describe the run.
+
+## Dependencies
+
+None.
+
+## Implementation Approach
+
+A separate trait rather than more methods on `Integrator`. `Integrator` already
+carries `const ADAPTIVE: bool` because not every method is adaptive, and a
+fixed-step method has no tolerances to set. Splitting them keeps that honest and
+leaves room for a fixed-step integrator that implements `Integrator` and not
+this.
+
+```rust
+/// An adaptive integrator's error tolerances.
+pub trait Tolerances<const D: usize>: Integrator<D> {
+    /// The absolute tolerance, per component.
+    fn absolute_tolerance(&self) -> SVector<f64, D>;
+
+    /// The relative tolerance.
+    fn relative_tolerance(&self) -> f64;
+
+    /// Sets the absolute tolerance for every component.
+    fn a_tol(self, a_tol: f64) -> Self;
+
+    /// Sets the absolute tolerance per component.
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self;
+
+    /// Sets the relative tolerance.
+    fn r_tol(self, r_tol: f64) -> Self;
+}
+```
+
+The setters keep the names they already have, so **no existing call site
+changes**. `Vern7::new().a_tol(1e-12).r_tol(1e-12)` continues to compile,
+resolving through the trait instead of the inherent impl, because `Integrator` is
+already in the prelude and `Tolerances` would be too.
+
+The getters need different names because `a_tol` is taken by the setter. Naming
+them `absolute_tolerance`/`relative_tolerance` avoids the collision and reads
+better than `a_tol_value`.
+
+### Steps
+
+1. Declare `Tolerances<D>` in `src/integrator/mod.rs`.
+2. For each of `Vern7`, `DormandPrince45`, `BS3`, `Rosenbrock23`: move the three
+   setter bodies from the inherent impl into an `impl Tolerances<D>`, and add the
+   two getters returning the stored fields. Keep `new()` inherent.
+3. Export `Tolerances` from the prelude.
+4. Delete `astrokit`'s local copy and import this one.
+
+### What not to do
+
+Do not add the tolerance methods to `Integrator` itself. That would oblige a
+fixed-step method to implement setters that mean nothing, and `const ADAPTIVE`
+exists precisely to admit such methods.
+
+Do not add `Default` or `Copy` as supertraits of `Integrator` while doing this,
+tempting as it is for downstream bounds. All four current integrators are both,
+but a multistep method holding Nordsieck history would reasonably be neither, and
+a downstream crate can bound on them explicitly where it needs them.
+
+## Testing Requirements
+
+- Every integrator round-trips: set a tolerance, read it back.
+- `a_tol` sets all components; `a_tol_full` sets them independently.
+- A generic function over `S: Tolerances<D>` compiles and changes step counts as
+  the tolerance changes, for each integrator.
+- Existing tests continue to pass unchanged, which is the point: the setter names
+  do not move.
+
+## Success Criteria
+
+- A downstream crate generic over `S: Tolerances<D>` can set and read tolerances
+  for every adaptive integrator in this crate.
+- No call site in this crate or its tests changes.
+- `astrokit` deletes its local `Tolerances`.
+
+**Effort**: Small

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

@@ -0,0 +1,178 @@
+# Feature: Runtime Integrator Selection
+
+**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
+
+Which integrator runs a problem is fixed at compile time, because `Integrator<D>`
+keeps its structural facts in associated constants:
+
+```rust
+pub trait Integrator<const D: usize> {
+    const ORDER: usize;
+    const STAGES: usize;
+    const ADAPTIVE: bool;
+    const DENSE: bool;
+    const MAIN_STAGES: usize = Self::STAGES;
+    const EXTRA_STAGES: usize = 0;
+    const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES;
+}
+```
+
+`Problem` is generic over `S`, and `Problem::solve` reads `S::ADAPTIVE` to choose
+its stepping branch and `S::EXTRA_STAGES` / `S::TOTAL_DENSE_STAGES` to decide
+whether lazy stages still need computing. A value cannot answer those questions —
+only a type can.
+
+So the obvious workaround does not work. An enum wrapping the four integrators
+cannot implement `Integrator<D>`, because it would have to give one
+`EXTRA_STAGES` for all of them, and Vern7 has six where Dormand-Prince has none.
+
+`Box<dyn Integrator<D>>` is not available either, and will not become available:
+`step<T: Scalar, P>` and `interpolate<T: Scalar>` are generic methods, so the
+trait is not object-safe. Erasing the scalar would defeat the genericity that
+makes autodiff work. The enum is the practical route, not `dyn`.
+
+## Why This Feature Matters
+
+This surfaced building `astrokit` on top of this crate. Three things a dependent
+wants and cannot have:
+
+1. **Accuracy presets that include a method.** `Accuracy::High` should mean "tight
+   tolerances *and* a high-order method", since the two belong together — a 7th
+   order method at 1e-6 is waste, and a 3rd order method at 1e-12 is hopeless. As
+   it stands a preset can only set tolerances, and the method has to be named
+   separately as a type parameter.
+
+2. **An integrator named in a configuration file.** `astrokit` is built for a
+   study where every run parameter must be declared in a committed file and every
+   propagation must record enough to recreate itself. `integrator: "Vern7"`
+   deserialized from that file cannot select anything, because the choice is a
+   type. The integrator ends up recorded in the provenance but not settable from
+   it, which is backwards.
+
+3. **A heterogeneous collection.** Benchmarking or comparing methods over the same
+   problem means either a macro over method types or duplicated code per method,
+   where a `Vec<Method>` would do.
+
+None of these is fatal. `astrokit` has settled on Vern7 as a type-level default,
+which is the right method for its accuracy range anyway. But the constraint is
+structural rather than deliberate, and it is worth fixing before more dependents
+grow around it.
+
+## Dependencies
+
+None, though it interacts with [Tolerances Trait](40-tolerances-trait.md): an
+enum wrapper needs to forward tolerance setters too, which is much easier once
+those are trait methods.
+
+## Implementation Approach
+
+### Step 1: constants become methods
+
+The structural facts move from associated constants to methods taking `&self`:
+
+```rust
+pub trait Integrator<const D: usize> {
+    fn order(&self) -> usize;
+    fn stages(&self) -> usize;
+    fn adaptive(&self) -> bool;
+    fn dense(&self) -> bool;
+
+    fn main_stages(&self) -> usize { self.stages() }
+    fn extra_stages(&self) -> usize { 0 }
+    fn total_dense_stages(&self) -> usize { self.main_stages() + self.extra_stages() }
+
+    // step, interpolate, compute_extra_stages unchanged
+}
+```
+
+**This costs nothing measurable.** Every use of these constants today is already a
+runtime comparison, not a compile-time size:
+
+- `if S::ADAPTIVE` in `Problem::solve` becomes `if self.integrator.adaptive()`,
+  which monomorphizes to the same constant and folds the same way for a concrete
+  integrator.
+- `borrowed.len() < 2 + S::TOTAL_DENSE_STAGES` compares against a `Vec` length.
+- `&dense[2..2 + S::MAIN_STAGES]` slices a `Vec`.
+
+None of them appears in type position, so nothing needs the values to be consts.
+For a concrete `S` the method calls inline to the same constants they are now.
+
+This is a breaking change to anyone implementing `Integrator` by hand — but a
+mechanical one, and this crate has no external implementors.
+
+### Step 2: an enum in this crate
+
+With step 1 done, an enum can dispatch:
+
+```rust
+/// An integrator chosen at runtime.
+#[non_exhaustive]
+#[derive(Debug, Clone, Copy)]
+pub enum Method<const D: usize> {
+    BS3(BS3<D>),
+    DormandPrince45(DormandPrince45<D>),
+    Vern7(Vern7<D>),
+    Rosenbrock23(Rosenbrock23<D>),
+}
+
+impl<const D: usize> Integrator<D> for Method<D> {
+    fn extra_stages(&self) -> usize {
+        match self {
+            Method::Vern7(m) => m.extra_stages(),
+            Method::BS3(m) => m.extra_stages(),
+            ..
+        }
+    }
+    // and so on, including step and interpolate
+}
+```
+
+Each method body is a `match` forwarding to the wrapped integrator. Verbose, and
+a good candidate for a macro over the variant list.
+
+Worth deriving `Serialize`/`Deserialize` on a companion `MethodName` enum so a
+method can be named in a configuration file, and adding
+`MethodName::instantiate(&self) -> Method<D>`. That is what makes point 2 above
+work.
+
+### Step 3: presets
+
+Once a method is a value, an accuracy preset can carry one, pairing a tolerance
+with a method whose order suits it. That belongs in whichever crate owns the
+notion of a preset — probably this one, given
+[Default Algorithm Selection](13-default-algorithm-selection.md) is already on
+the roadmap and needs the same machinery.
+
+### The cost of not doing step 1
+
+If the constants must stay, the fallback is a parallel object-ish layer: a second
+trait with `&self` methods, blanket-implemented for everything that implements
+`Integrator`, and `Problem` generic over that instead. It works, it duplicates the
+trait surface, and every implementor carries both. Step 1 is cheaper and cleaner.
+
+## Testing Requirements
+
+- Every existing test passes after step 1, with method calls substituted for
+  constant reads.
+- A benchmark confirms no regression from the constants becoming methods — the
+  claim is that they inline, and that should be measured rather than assumed, on
+  the existing `simple_1d` and `orbit` benches.
+- `Method` gives bit-identical results to the integrator it wraps, for each
+  variant, on a problem with dense output interpolation (which exercises the
+  lazy-stage path where `extra_stages` matters).
+- A method name round-trips through serde and instantiates the right variant.
+
+## Success Criteria
+
+- `Problem` can be constructed with an integrator chosen at runtime.
+- No measurable performance change on the existing criterion benches.
+- A dependent crate can offer an accuracy preset that selects both a tolerance and
+  a method, and can read a method name out of a configuration file.
+
+**Effort**: Medium — step 1 is small and mechanical, step 2 is repetitive

roadmap/features/42-callback-effects.md

@@ -0,0 +1,231 @@
+# Feature: Callback Effects That Can Change the State
+
+**Status**: completed (2026-10-01)
+
+## Overview
+
+Continuous callbacks located their crossings correctly and then threw away what
+the effect did about them. In `Problem::solve`:
+
+```rust
+callback.effect(&mut self.ode);   // writes to ode.y
+}
+}
+self.ode.y = new_y;               // and here it goes
+self.ode.t += curr_step;
+```
+
+The effect ran on the pre-step state, and the two lines after the callback loop
+overwrote `ode.y` with the state the integrator had produced. Anything the effect
+wrote to the state was lost. `stop` worked only because it writes `t_end`
+instead.
+
+Reproduced with `dy/dt = 1` from 0 to 10, `y0 = 0`, an event of `5 - y[0]` and an
+effect of `ode.y = Vector1::new(100.0)`: the solution ended at `y = 10`, exactly
+as it does with an empty effect, where it should end near 105.
+
+So the machinery was in place for the easy half of event handling — detect a
+condition and stop — and could not do the other half: detect a condition, change
+the state, and carry on.
+
+## Why This Feature Matters
+
+A rotating rigid body parameterized by Modified Rodrigues Parameters has a
+singularity at a full turn. It is avoided by swapping the parameter vector for
+the shadow set `-sigma / |sigma|^2` — the same attitude, different parameters —
+whenever `|sigma|` exceeds 1. The event is `1 - |sigma|`, the effect is that
+substitution, and the integration has to continue from the substituted state.
+Any arc longer than a half turn crosses it, so this is the common case rather
+than an edge case.
+
+The same shape of problem covers a body entering or leaving shadow, a vehicle
+whose mass or area changes partway through a run, and a bouncing ball. All of
+them are a state discontinuity located by an event, which is the one thing the
+callbacks could not express.
+
+## What Was Done
+
+### The effect runs after the crossing is committed
+
+The event branch now commits the shortened step before calling the effect:
+
+1. the step is shortened to land on the crossing, and the state there is
+   recorded;
+2. `ode.t` and `ode.y` are set to that crossing;
+3. the effect is called, with the ODE at the crossing;
+4. whatever the effect leaves in `ode.y` is what the integration continues from.
+
+Detection still reads the state at the start of the step — an event function is
+asked about both ends — so nothing can be committed until every callback has been
+asked. That is why this is a reordering of the branch rather than a swap of two
+lines.
+
+### `effect(&mut self)`
+
+`Callback::effect` takes `&mut self`, so a callback can record that it has fired,
+count its firings, or keep anything else it needs across them. `Problem` holds
+`Vec<&'a mut dyn Callback<T, D, P>>` and `with_callback` takes `&'a mut dyn
+Callback`, which is the breaking part of the change: a registered callback has to
+be declared `mut` and passed as `&mut`.
+
+`FnCallback` still holds its two closures by shared reference and so still has no
+state of its own; a callback that needs state implements `Callback` on a struct,
+which is what [Discrete Callbacks](05-discrete-callbacks.md) will do with its
+`has_fired` field.
+
+`Problem` lost its `#[derive(Clone)]`, which a `&mut` field cannot support. It
+was never used.
+
+### A terminal effect ends the solution at the crossing
+
+`ODE::finished` existed, was never set and was never read. It now means what its
+name says: `ODE::terminate` sets it, the solve loop checks it, and `stop` is
+`ode.terminate()`.
+
+Behaviorally this is close to what `stop` already did. Writing `t_end = ode.t`
+also ended the loop, and because the step had already been shortened to the
+crossing, the last recorded point was the crossing — the solution did not in fact
+run a step past the event. What the old `stop` did do was mutate the problem's own
+end time, so a solution could no longer report the `t_end` it was asked for. Now
+`t_end` is left alone and `finished` says why the solution is short.
+
+### Both sides of a jump are recorded
+
+Where an effect changes the state, the trajectory genuinely has two values at
+that time, and the solution holds both: the same `t` appears twice in `times`,
+with the state at the crossing and the state the effect left. The `dense` entry
+between them spans zero time and is never interpolated into — `interpolate` at
+exactly that time hits the search's exact-match branch, and walks to the last
+entry at that time, so it reports the state the effect left. That is the value the
+trajectory carries forward, and it is the one a caller asking "where is the body
+at this instant" wants; the value it arrived with is the entry before it in
+`states`.
+
+This is what DifferentialEquations.jl does, where the two are `u_left` and
+`u_right`.
+
+### The controller is reset across a jump
+
+`Controller` grew a required `reset` method, which discards the error history and
+leaves the next step guess alone. `solve` calls it when an effect changes the
+state. The history is a claim about how the error behaves along a smooth
+solution, and across a discontinuity it is not about anything.
+
+Keyed on the state changing rather than on the effect running at all, so that a
+callback which only observes does not throw away a good history every time it
+fires. An effect that changes only the parameters is a discontinuity the solver
+cannot see, and does not reset the controller.
+
+`PIDController::reset` already existed as an inherent method and became the trait
+implementation; `PIController::reset` restores `factor_old` to its initial value.
+
+### Two events in one step
+
+Defined and tested: the earliest crossing in the step fires, the step ends there,
+and any later event is still ahead and is found again on the following step. A
+tie goes to registration order.
+
+The old behavior was to walk the callbacks in registration order, each shortening
+what the last one left. That could apply a later event's effect and then truncate
+the step to before it, so the effect was applied at a time the solution never
+reached. One pass to find the earliest crossing costs an extra root find only
+when more than one callback crosses in the same step.
+
+### Crossings at exactly zero, and not firing twice at one
+
+The test was `before * after < 0.0`, which misses a crossing that starts or ends
+exactly on zero. Now:
+
+- `before == 0.0` is a crossing at the step's own start, located at `h = 0`;
+- `after == 0.0` is a crossing at the step's end;
+- opposite signs are located with Brent as before.
+
+The first of those needs a guard, because an effect is free to leave its event
+function sitting on the zero it fired at. The shadow set swap does exactly that:
+`|sigma|` is 1 on both sides of it, so `1 - |sigma|` is still zero after the
+effect. Without a guard the callback is re-detected at the start of the next
+step, located at `h = 0`, fired again, forever.
+
+The guard is that the step beginning at a firing searches from
+`EVENT_PROBE_FRACTION` of itself — a millionth — rather than from its start, for
+that callback alone. Both the `before` sign and Brent's bracket start there, so
+the event just handled is out of reach while a second crossing later in the same
+step is still found. `solve` tracks the time each callback last fired to know
+when to do this.
+
+The first attempt at this guard was to skip a callback for the whole step that
+begins at its firing, which is simpler and wrong: a bouncing ball bounced once
+and then fell through the floor, because the step after the bounce contained the
+next landing and the callback was not looking. That is in the test suite now.
+
+A fraction of a step is not a satisfying thing to have to choose, and this is the
+least bad option found rather than a good one. Alternatives considered: an
+absolute tolerance on the event value, which is unit-dependent and no better
+founded; taking the sign from the step before the firing, which still leaves
+Brent bracketing a root at `h = 0`; and refusing a located crossing below some
+`h`, which is the same arbitrary fraction one step further down. If something
+cleaner turns up, this is the place to change it.
+
+## Known Limitations
+
+- An even number of crossings inside one step is still missed: the two ends of
+  the step agree in sign and nothing looks in between. Documented on
+  `Problem::solve`.
+- A crossing inside the first millionth of the step that begins at a firing is
+  missed, which is the cost of the no-refiring guard. It takes an event that
+  recurs almost immediately after one was handled.
+- The crossing time is located in real arithmetic and carries no derivative.
+  Deliberate: differentiating through Brent's iterations gives the derivative of
+  the iterative procedure rather than the sensitivity of the event, which is
+  wrong without being loud about it. The state at the crossing comes from one
+  ordinary step of the integrator, so it does carry derivatives.
+- An effect that changes the parameters rather than the state does not reset the
+  controller, and does not record two points.
+
+## Testing
+
+In `src/problem.rs`:
+
+- the reproduction above, ending near 105 rather than 10;
+- both sides of the jump present in `times`/`states`, interpolation either side of
+  it following the branch it is on, and interpolation at the jump itself giving
+  the state the effect left;
+- `stop` ending the solution at the crossing, with nothing past it, `t_end`
+  untouched and `finished` set;
+- a callback counting its own firings and acting only on the first, which needs
+  the `&mut self`;
+- two events in one step firing earliest-first, with the registration order
+  reversed to show it does not matter;
+- an effect that lands back on its own event firing once there rather than
+  looping, and firing again when the event genuinely recurs;
+- a bouncing ball, bouncing at the three analytic bounce times inside the
+  interval and interpolating correctly inside each arc;
+- an event function already zero at the start of the integration firing there;
+- the controller reset asked for once across a jump and not at all for an effect
+  that changes nothing, measured with a controller that counts the resets;
+- a jump under dual numbers, where the derivative through the effect is exact
+  (in `src/lib.rs`, with the autodiff tests);
+- the Modified Rodrigues Parameter swap itself, against the closed form for a
+  constant angular velocity: one swap, at exactly half a turn, `|sigma| <= 1`
+  throughout, and the final attitude equal to `-cot(theta / 4)` about the
+  rotation axis.
+
+## Success Criteria
+
+- An effect that writes to `ode.y` changes what is integrated from. ✅
+- An effect can record something about itself. ✅
+- A terminated solution ends at the crossing. ✅
+- A jump is recorded on both sides. ✅
+- The step size controller forgets its history across a jump. ✅
+- The order of two events in one step is defined and tested. ✅
+- No measurable change on the criterion benches. ✅ Measured before and after on
+  the same machine in one sitting: `bench_orbit` -1.02%, `bench_simple_1d`
+  +0.36%, `bench_interpolation_1d` -0.67%, all three within criterion's noise
+  threshold. Which is the right reading of them: re-running the *identical*
+  binary moved `bench_simple_1d` by -0.80% and `bench_interpolation_1d` by
+  +1.47%, so run-to-run drift here is about a percent either way and is larger
+  than the difference being looked for. No bench registers a callback, so what
+  they measure is the no-callback path, which gains a loop over an empty slice
+  and a `finished` check per step.
+
+**Effort**: Medium

src/callback.rs

@@ -1,31 +1,96 @@
+//! Conditions that interrupt an integration, and what they do when they fire.
+
 use nalgebra::SVector;
 
 use crate::ode::ODE;
 use crate::scalar::Scalar;
 
-/// A callback function you can utilize to change your model mid-run without stopping.
+/// A condition the solver watches for, and the change it makes when it is met.
 ///
 /// The solver evaluates [`event`](Callback::event) at each end of every accepted
-/// step and looks for a sign change. When it finds one it locates the crossing
-/// and applies [`effect`](Callback::effect) to the ODE there.
+/// step and looks for a sign change. When it finds one it shortens the step to
+/// land on the crossing, records the state there, and then calls
+/// [`effect`](Callback::effect) on the ODE. Whatever the effect leaves in
+/// `ode.y` is the state the integration carries on from.
+///
+/// An effect may:
+///
+/// - change the state, which records a jump: the solution holds both the value
+///   at the crossing and the value the effect left, at the same time (see
+///   [`Solution`](crate::problem::Solution));
+/// - change the parameters, which changes the right-hand side from the crossing
+///   onward;
+/// - end the solution at the crossing, by calling [`stop`] or
+///   [`ODE::terminate`];
+/// - record something about itself, since it takes `&mut self`.
+///
+/// # Example
+///
+/// A fire-once callback, which needs the `&mut self` on the effect to remember
+/// that it has run:
+///
+/// ```
+/// use nalgebra::Vector1;
+/// use ordinary_diffeq::prelude::*;
+///
+/// struct Kick {
+///     fired: usize,
+/// }
+///
+/// impl Callback<f64, 1, ()> for Kick {
+///     fn event(&self, _t: f64, y: &Vector1<f64>, _p: &()) -> f64 {
+///         5.0 - y[0]
+///     }
+///
+///     fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) {
+///         self.fired += 1;
+///         ode.y = Vector1::new(100.0);
+///     }
+/// }
+///
+/// fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+///     Vector1::new(1.0)
+/// }
+///
+/// let ode = ODE::new(&derivative, 0.0, 10.0, Vector1::new(0.0), ());
+/// let mut kick = Kick { fired: 0 };
+/// let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+///     .with_callback(&mut kick);
+/// let solution = problem.solve();
+///
+/// // The state jumped to 100 at t = 5 and kept integrating from there.
+/// assert!((solution.states.last().unwrap()[0] - 105.0).abs() < 1e-9);
+/// ```
+pub trait Callback<T: Scalar, const D: usize, P> {
+    /// The function whose zero crossings mark the event.
+    ///
+    /// Called repeatedly while the crossing is being located, so it must not
+    /// depend on how many times it has been called.
+    fn event(&self, t: T, y: &SVector<T, D>, params: &P) -> T;
+
+    /// What to do to the ODE once the crossing has been located.
+    ///
+    /// `ode.t` and `ode.y` are the time and state at the crossing when this is
+    /// called.
+    fn effect(&mut self, ode: &mut ODE<'_, T, D, P>);
+}
+
+/// A [`Callback`] built from a pair of closures.
+///
+/// The closures are held by shared reference, so neither can carry state of its
+/// own. A callback that has to remember something — that it has already fired,
+/// how many times, what the state was — implements [`Callback`] on a struct
+/// instead; the example on that trait shows it.
+///
 /// ```
 /// # use ordinary_diffeq::prelude::*;
 /// # use nalgebra::Vector2;
 /// # type Params = ();
-/// let too_high = FnCallback {
+/// let mut too_high = FnCallback {
 ///     event: &|_t: f64, y: &Vector2<f64>, _p: &Params| 10.0 - y[0],
 ///     effect: &stop,
 /// };
 /// ```
-pub trait Callback<T: Scalar, const D: usize, P> {
-    /// The function whose zero crossings mark the event.
-    fn event(&self, t: T, y: &SVector<T, D>, params: &P) -> T;
-
-    /// What to do to the ODE once the crossing has been located.
-    fn effect(&self, ode: &mut ODE<'_, T, D, P>);
-}
-
-/// A [`Callback`] built from a pair of closures.
 #[derive(Clone, Copy)]
 pub struct FnCallback<'a, T: Scalar, const D: usize, P> {
     /// The function to check for zero crossings
@@ -42,14 +107,14 @@ impl<T: Scalar, const D: usize, P> Callback<T, D, P> for FnCallback<'_, T, D, P>
     }
 
     #[inline]
-    fn effect(&self, ode: &mut ODE<'_, T, D, P>) {
+    fn effect(&mut self, ode: &mut ODE<'_, T, D, P>) {
         (self.effect)(ode)
     }
 }
 
-/// A convenience function for stopping the integration
+/// An effect that ends the solution at the crossing.
 pub fn stop<T: Scalar, const D: usize, P>(ode: &mut ODE<'_, T, D, P>) {
-    ode.t_end = ode.t;
+    ode.terminate();
 }
 
 #[cfg(test)]

src/controller.rs

@@ -34,6 +34,14 @@ pub trait Controller<T: Scalar> {
     fn next_step_guess(&self) -> TryStep<T>;
 
     fn set_next_step_guess(&mut self, guess: TryStep<T>);
+
+    /// Discards the error history, leaving the next step guess alone.
+    ///
+    /// The history is an estimate of how the error behaves along a smooth
+    /// solution, so it says nothing across a discontinuity.
+    /// [`Problem::solve`](crate::problem::Problem::solve) calls this when a
+    /// callback effect changes the state.
+    fn reset(&mut self);
 }
 
 #[derive(Debug, Clone, Copy)]
@@ -68,9 +76,7 @@ impl<T: Scalar> Controller<T> for PIController<T> {
             TryStep::Accepted(prev_step, h)
         } else {
             // Reject the stepsize and propose a smaller one for the current step
-            TryStep::NotYetAccepted(
-                prev_step / self.factor_c1.min(factor_11 / self.safety_factor),
-            )
+            TryStep::NotYetAccepted(prev_step / self.factor_c1.min(factor_11 / self.safety_factor))
         }
     }
 
@@ -82,6 +88,9 @@ impl<T: Scalar> Controller<T> for PIController<T> {
         self.next_step_guess = guess;
     }
 
+    fn reset(&mut self) {
+        self.factor_old = 1.0e-4;
+    }
 }
 
 impl<T: Scalar> PIController<T> {
@@ -113,7 +122,6 @@ impl<T: Scalar> Default for PIController<T> {
     }
 }
 
-
 /// PID (Proportional-Integral-Derivative) step size controller
 ///
 /// The PID controller is an advanced adaptive time-stepping controller that provides
@@ -160,20 +168,20 @@ impl<T: Scalar> Default for PIController<T> {
 #[derive(Debug, Clone, Copy)]
 pub struct PIDController<T: Scalar = f64> {
     // PID Coefficients
-    pub p_term: f64,  // Proportional: reaction to current error
-    pub i_term: f64,  // Integral: reaction to error history
-    pub d_term: f64,  // Derivative: reaction to error rate of change
+    pub p_term: f64, // Proportional: reaction to current error
+    pub i_term: f64, // Integral: reaction to error history
+    pub d_term: f64, // Derivative: reaction to error rate of change
 
     // Constraints
-    pub factor_c1: f64,      // 1/min_factor (maximum step decrease)
-    pub factor_c2: f64,      // 1/max_factor (maximum step increase)
-    pub h_max: f64,          // Maximum allowed step size
-    pub safety_factor: f64,  // Safety factor (typically 0.9)
+    pub factor_c1: f64,     // 1/min_factor (maximum step decrease)
+    pub factor_c2: f64,     // 1/max_factor (maximum step increase)
+    pub h_max: f64,         // Maximum allowed step size
+    pub safety_factor: f64, // Safety factor (typically 0.9)
 
     // Error history for PID control
-    pub err_old: f64,     // ε_{n-1}: previous step error
-    pub err_older: f64,   // ε_{n-2}: error two steps ago
-    pub h_old: f64,       // h_{n-1}: previous step size
+    pub err_old: f64,   // ε_{n-1}: previous step error
+    pub err_older: f64, // ε_{n-2}: error two steps ago
+    pub h_old: f64,     // h_{n-1}: previous step size
 
     // Next step guess
     pub next_step_guess: TryStep<T>,
@@ -188,15 +196,15 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
         let factor = if self.err_old <= 0.0 {
             // First step: use only proportional control
             let factor_11 = err.powf(self.p_term);
-            self.factor_c2.max(
-                self.factor_c1.min(factor_11 / self.safety_factor)
-            )
+            self.factor_c2
+                .max(self.factor_c1.min(factor_11 / self.safety_factor))
         } else if self.err_older <= 0.0 {
             // Second step: use PI control (proportional + integral)
             let factor_11 = err.powf(self.p_term);
             let factor_12 = self.err_old.powf(-self.i_term);
             self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 / self.safety_factor)
+                self.factor_c1
+                    .min(factor_11 * factor_12 / self.safety_factor),
             )
         } else {
             // Full PID control (proportional + integral + derivative)
@@ -209,7 +217,8 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
                 1.0
             };
             self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 * factor_13 / self.safety_factor)
+                self.factor_c1
+                    .min(factor_11 * factor_12 * factor_13 / self.safety_factor),
             )
         };
 
@@ -219,7 +228,7 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
 
             // Update error history for next step
             self.err_older = self.err_old;
-            self.err_old = err.max(1.0e-4);  // Prevent very small values
+            self.err_old = err.max(1.0e-4); // Prevent very small values
             self.h_old = h.re();
 
             // Apply maximum step size limit
@@ -247,6 +256,11 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
         self.next_step_guess = guess;
     }
 
+    fn reset(&mut self) {
+        self.err_old = 0.0;
+        self.err_older = 0.0;
+        self.h_old = 0.0;
+    }
 }
 
 impl<T: Scalar> PIDController<T> {
@@ -281,9 +295,9 @@ impl<T: Scalar> PIDController<T> {
             factor_c2: 1.0 / max_factor,
             h_max: h_max.abs(),
             safety_factor,
-            err_old: 0.0,      // No history initially
-            err_older: 0.0,    // No history initially
-            h_old: 0.0,        // No history initially
+            err_old: 0.0,   // No history initially
+            err_older: 0.0, // No history initially
+            h_old: 0.0,     // No history initially
             next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)),
         }
     }
@@ -301,25 +315,16 @@ impl<T: Scalar> PIDController<T> {
     pub fn for_order(order: usize) -> Self {
         let k_plus_1 = (order + 1) as f64;
         Self::new(
-            0.49 / k_plus_1,  // p_term: proportional
-            0.34 / k_plus_1,  // i_term: integral
-            0.10 / k_plus_1,  // d_term: derivative
-            10.0,              // max_factor
-            0.2,               // min_factor
-            100000.0,          // h_max
-            0.9,               // safety_factor
-            1e-4,              // initial_h
+            0.49 / k_plus_1, // p_term: proportional
+            0.34 / k_plus_1, // i_term: integral
+            0.10 / k_plus_1, // d_term: derivative
+            10.0,            // max_factor
+            0.2,             // min_factor
+            100000.0,        // h_max
+            0.9,             // safety_factor
+            1e-4,            // initial_h
         )
     }
-
-    /// Reset the controller's error history
-    ///
-    /// Useful when switching algorithms or restarting integration
-    pub fn reset(&mut self) {
-        self.err_old = 0.0;
-        self.err_older = 0.0;
-        self.h_old = 0.0;
-    }
 }
 
 impl<T: Scalar> Default for PIDController<T> {
@@ -335,19 +340,18 @@ impl<T: Scalar> Default for PIDController<T> {
     /// For method-specific tuning, use `PIDController::for_order(order)` instead.
     fn default() -> Self {
         Self::new(
-            0.07,      // p_term (proportional)
-            0.04,      // i_term (integral)
-            0.01,      // d_term (derivative)
-            10.0,      // max_factor
-            0.2,       // min_factor
-            100000.0,  // h_max
-            0.9,       // safety_factor
-            1e-4,      // initial_h
+            0.07,     // p_term (proportional)
+            0.04,     // i_term (integral)
+            0.01,     // d_term (derivative)
+            10.0,     // max_factor
+            0.2,      // min_factor
+            100000.0, // h_max
+            0.9,      // safety_factor
+            1e-4,     // initial_h
         )
     }
 }
 
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -399,7 +403,7 @@ mod tests {
         // Default uses conservative PID coefficients
         assert_eq!(controller.p_term, 0.07);
         assert_eq!(controller.i_term, 0.04);
-        assert_eq!(controller.d_term, 0.01);  // True PID with derivative term
+        assert_eq!(controller.d_term, 0.01); // True PID with derivative term
     }
 
     #[test]
@@ -436,7 +440,8 @@ mod tests {
 
         for &err in &errors {
             let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
+            let mut pid_result =
+                <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
 
             if pi_result.is_accepted() {
                 pi_steps.push(pi_result.extract());
@@ -548,7 +553,8 @@ mod tests {
 
         for &err in &errors {
             let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
+            let mut pid_result =
+                <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
 
             if pi_result.is_accepted() {
                 let pi_h_next = pi_result.reset().unwrap().extract();
@@ -567,11 +573,8 @@ mod tests {
 
         // Compute standard deviation of log step size ratios
         let pi_mean: f64 = pi_ratios.iter().sum::<f64>() / pi_ratios.len() as f64;
-        let pi_variance: f64 = pi_ratios
-            .iter()
-            .map(|r| (r - pi_mean).powi(2))
-            .sum::<f64>()
-            / pi_ratios.len() as f64;
+        let pi_variance: f64 =
+            pi_ratios.iter().map(|r| (r - pi_mean).powi(2)).sum::<f64>() / pi_ratios.len() as f64;
         let pi_std = pi_variance.sqrt();
 
         let pid_mean: f64 = pid_ratios.iter().sum::<f64>() / pid_ratios.len() as f64;

src/integrator/bs3.rs

@@ -1,9 +1,12 @@
 use nalgebra::SVector;
 
 use super::super::ode::ODE;
-use super::{Integrator, StepResult};
+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];
@@ -61,7 +64,7 @@ pub struct BS3<const D: usize> {
 impl<const D: usize> Default for BS3<D>
 where
     BS3<D>: Integrator<D>,
- {
+{
     fn default() -> Self {
         Self::new()
     }
@@ -80,21 +83,31 @@ where
             r_tol: 1e-8,
         }
     }
+}
 
-    /// Set absolute tolerance (same value for all components)
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for BS3<D>
+where
+    BS3<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
         self.a_tol = SVector::<f64, D>::from_element(a_tol);
         self
     }
 
-    /// Set absolute tolerance (different value per component)
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
         self.a_tol = a_tol;
         self
     }
 
-    /// Set relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
+    fn r_tol(mut self, r_tol: f64) -> Self {
         self.r_tol = r_tol;
         self
     }
@@ -107,6 +120,10 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     // Row 2: [1/2]
     // Row 3: [0, 3/4]
     // Row 4: [2/9, 1/3, 4/9]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
     const A: &'a [f64] = &[
         1.0 / 2.0,           // a[1,0]
         0.0,                 // a[2,0]
@@ -117,6 +134,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Solution weights (3rd order)
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const B: &'a [f64] = &[
         2.0 / 9.0,           // b[0]
         1.0 / 3.0,           // b[1]
@@ -125,6 +144,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Error estimate weights (difference between 3rd and 2nd order)
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const B_ERROR: &'a [f64] = &[
         2.0 / 9.0 - 7.0 / 24.0,      // b[0] - b*[0]
         1.0 / 3.0 - 1.0 / 4.0,       // b[1] - b*[1]
@@ -133,6 +154,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Stage times
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const C: &'a [f64] = &[
         0.0,                 // c[0]
         1.0 / 2.0,           // c[1]
@@ -145,14 +168,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);
@@ -167,7 +201,9 @@ where
 
         // Stage 4: k4 = f(t + c[3]*h, y + h*(a[3,0]*k1 + a[3,1]*k2 + a[3,2]*k3))
         let y4 = ode.y
-            + (k[0] * T::lift(Self::A[3]) + k[1] * T::lift(Self::A[4]) + k[2] * T::lift(Self::A[5]))
+            + (k[0] * T::lift(Self::A[3])
+                + k[1] * T::lift(Self::A[4])
+                + k[2] * T::lift(Self::A[5]))
                 * h;
         k[3] = (ode.f)(ode.t + h * Self::C[3], y4, &ode.params);
 
@@ -201,10 +237,10 @@ where
         // BS3 uses standard cubic Hermite interpolation with derivatives at endpoints
         // Store: y0, y1, f0=k[0], f1=k[3] (FSAL)
         let dense_coeffs = vec![
-            ode.y,           // y0 at start of step
-            next_y,          // y1 at end of step
-            k[0],            // f(t0, y0) - derivative at start
-            k[3],            // f(t1, y1) - derivative at end (FSAL)
+            ode.y,  // y0 at start of step
+            next_y, // y1 at end of step
+            k[0],   // f(t0, y0) - derivative at start
+            k[3],   // f(t1, y1) - derivative at end (FSAL)
         ];
 
         StepResult {
@@ -261,13 +297,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]
@@ -281,7 +316,7 @@ mod tests {
         let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
 
         let bs3 = BS3::new();
-        let h = 0.001;  // Smaller step size for tighter tolerances
+        let h = 0.001; // Smaller step size for tighter tolerances
 
         let __step = bs3.step(&ode, h);
         let (y_next, err, dense) = (__step.y, __step.err, __step.dense);
@@ -311,7 +346,7 @@ mod tests {
         let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
 
         let bs3 = BS3::new();
-        let h = 0.001;  // Smaller step size
+        let h = 0.001; // Smaller step size
 
         let __step = bs3.step(&ode, h);
         let (_y_next, _err, dense) = (__step.y, __step.err, __step.dense);

src/integrator/dormand_prince.rs

@@ -1,9 +1,13 @@
 use nalgebra::SVector;
 
-use super::{Integrator, StepResult};
+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];
@@ -73,25 +77,35 @@ where
             r_tol: 1e-8,
         }
     }
+}
 
-    /// The same absolute tolerance on every component.
-    pub fn a_tol(self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for DormandPrince45<D>
+where
+    DormandPrince45<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(self, a_tol: f64) -> Self {
         Self {
             a_tol: SVector::<f64, D>::from_element(a_tol),
             r_tol: self.r_tol,
         }
     }
 
-    /// A per-component absolute tolerance.
-    pub fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
         Self {
             a_tol,
             r_tol: self.r_tol,
         }
     }
 
-    /// The same relative tolerance on every component
-    pub fn r_tol(self, r_tol: f64) -> Self {
+    fn r_tol(self, r_tol: f64) -> Self {
         Self {
             a_tol: self.a_tol,
             r_tol,
@@ -164,39 +178,54 @@ 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) {
                 y_term += item * T::lift(Self::A[(i * (i - 1)) / 2 + j]);
             }
-            k[i] = (ode.f)(ode.t + T::lift(Self::C[i]) * h, ode.y + y_term * h, &ode.params);
+            k[i] = (ode.f)(
+                ode.t + T::lift(Self::C[i]) * h,
+                ode.y + y_term * h,
+                &ode.params,
+            );
 
             // 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>;
@@ -88,6 +110,46 @@ pub trait Integrator<const D: usize> {
     }
 }
 
+/// An adaptive integrator's error tolerances.
+///
+/// Split out from [`Integrator`] rather than folded into it, because a
+/// fixed-step method has no tolerances to set and `ADAPTIVE` exists to admit
+/// such methods. Every adaptive integrator in this crate implements it, so a
+/// function generic over `S: Tolerances<D>` can set and read the tolerances of
+/// whichever one it is handed.
+///
+/// # Example
+///
+/// ```rust
+/// use ordinary_diffeq::prelude::*;
+/// use nalgebra::SVector;
+///
+/// /// Tightens whatever integrator it is given.
+/// fn high_accuracy<const D: usize, S: Tolerances<D>>(integrator: S) -> S {
+///     integrator.a_tol(1e-12).r_tol(1e-12)
+/// }
+///
+/// let vern7 = high_accuracy(Vern7::<3>::new());
+/// assert_eq!(vern7.relative_tolerance(), 1e-12);
+/// assert_eq!(vern7.absolute_tolerance(), SVector::<f64, 3>::from_element(1e-12));
+/// ```
+pub trait Tolerances<const D: usize>: Integrator<D> {
+    /// The absolute tolerance, per component.
+    fn absolute_tolerance(&self) -> SVector<f64, D>;
+
+    /// The relative tolerance.
+    fn relative_tolerance(&self) -> f64;
+
+    /// Sets the absolute tolerance for every component.
+    fn a_tol(self, a_tol: f64) -> Self;
+
+    /// Sets the absolute tolerance per component.
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self;
+
+    /// Sets the relative tolerance.
+    fn r_tol(self, r_tol: f64) -> Self;
+}
+
 #[cfg(test)]
 mod tests {
     use super::dormand_prince::*;
@@ -119,3 +181,101 @@ mod tests {
         }
     }
 }
+
+#[cfg(test)]
+mod tolerance_tests {
+    use super::bs3::BS3;
+    use super::dormand_prince::DormandPrince45;
+    use super::rosenbrock::Rosenbrock23;
+    use super::vern7::Vern7;
+    use super::Tolerances;
+    use crate::controller::PIController;
+    use crate::problem::Problem;
+    use nalgebra::{SVector, Vector1};
+
+    /// Sets both tolerances and reads them back, for one integrator.
+    fn round_trip<S: Tolerances<3>>(integrator: S) {
+        let set = integrator.a_tol(1e-11).r_tol(1e-9);
+        assert_eq!(set.absolute_tolerance(), SVector::<f64, 3>::repeat(1e-11));
+        assert_eq!(set.relative_tolerance(), 1e-9);
+    }
+
+    /// `a_tol` sets every component; `a_tol_full` sets them independently.
+    fn per_component<S: Tolerances<3>>(integrator: S) {
+        let r_tol = integrator.relative_tolerance();
+        let uniform = integrator.a_tol(1e-7);
+        assert_eq!(
+            uniform.absolute_tolerance(),
+            SVector::<f64, 3>::repeat(1e-7)
+        );
+
+        let per = SVector::<f64, 3>::new(1e-6, 1e-9, 1e-12);
+        let full = uniform.a_tol_full(per);
+        assert_eq!(full.absolute_tolerance(), per);
+        // The relative tolerance is left alone by either setter.
+        assert_eq!(full.relative_tolerance(), r_tol);
+    }
+
+    #[test]
+    fn test_tolerances_round_trip() {
+        round_trip(BS3::<3>::new());
+        round_trip(DormandPrince45::<3>::new());
+        round_trip(Vern7::<3>::new());
+        round_trip(Rosenbrock23::<3>::new());
+    }
+
+    #[test]
+    fn test_tolerances_per_component() {
+        per_component(BS3::<3>::new());
+        per_component(DormandPrince45::<3>::new());
+        per_component(Vern7::<3>::new());
+        per_component(Rosenbrock23::<3>::new());
+    }
+
+    /// Solves `y' = -y` with whatever integrator it is handed, and reports how
+    /// many steps it took. Generic over `S`, which is the point of the trait.
+    fn steps_taken<S: Tolerances<1> + Copy>(integrator: S, a_tol: f64, r_tol: f64) -> usize {
+        fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+            -y
+        }
+
+        let ode = crate::ode::ODE::new(&derivative, 0.0, 10.0, Vector1::new(1.0), ());
+        let mut problem = Problem::new(
+            ode,
+            integrator.a_tol(a_tol).r_tol(r_tol),
+            PIController::default(),
+        );
+        problem.solve().times.len()
+    }
+
+    #[test]
+    fn test_generic_tightening_costs_steps() {
+        for (name, loose, tight) in [
+            (
+                "BS3",
+                steps_taken(BS3::<1>::new(), 1e-4, 1e-4),
+                steps_taken(BS3::<1>::new(), 1e-9, 1e-9),
+            ),
+            (
+                "DormandPrince45",
+                steps_taken(DormandPrince45::<1>::new(), 1e-4, 1e-4),
+                steps_taken(DormandPrince45::<1>::new(), 1e-9, 1e-9),
+            ),
+            (
+                "Vern7",
+                steps_taken(Vern7::<1>::new(), 1e-4, 1e-4),
+                steps_taken(Vern7::<1>::new(), 1e-9, 1e-9),
+            ),
+            (
+                "Rosenbrock23",
+                steps_taken(Rosenbrock23::<1>::new(), 1e-4, 1e-4),
+                steps_taken(Rosenbrock23::<1>::new(), 1e-9, 1e-9),
+            ),
+        ] {
+            assert!(
+                tight > loose,
+                "{name}: tightening the tolerance should cost steps, got {loose} loose and {tight} tight"
+            );
+        }
+    }
+}

src/integrator/rosenbrock.rs

@@ -1,7 +1,7 @@
 use nalgebra::{SMatrix, SVector};
 
 use super::super::ode::{ProblemFunction, ODE};
-use super::{Integrator, StepResult};
+use super::{Integrator, StepResult, Tolerances};
 use crate::scalar::Scalar;
 
 /// Strategy for when to update the Jacobian matrix
@@ -28,10 +28,11 @@ impl Default for JacobianUpdateStrategy {
 
 /// Compute the Jacobian matrix ∂f/∂y using forward finite differences
 ///
-/// For a system y' = f(t, y), computes the D×D Jacobian matrix J where J[i,j] = ∂f_i/∂y_j
+/// For a system `y' = f(t, y)`, computes the D×D Jacobian matrix `J` where
+/// `J[i,j] = ∂f_i/∂y_j`
 ///
-/// Uses forward differences: J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε
-/// where ε = √(machine_epsilon) * max(|y[j]|, 1.0)
+/// Uses forward differences: `J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε`
+/// where `ε = √(machine_epsilon) * max(|y[j]|, 1.0)`
 pub fn compute_jacobian<T: Scalar, const D: usize, P>(
     f: ProblemFunction<'_, T, D, P>,
     t: T,
@@ -143,8 +144,8 @@ pub struct Rosenbrock23<const D: usize> {
     c32: f64,
     /// Coefficient d = 1/(2+√2) ≈ 0.29289321881345254 (this is gamma!)
     d: f64,
-    /// Absolute tolerance for error estimation
-    a_tol: f64,
+    /// Absolute tolerance for error estimation, per component
+    a_tol: SVector<f64, D>,
     /// Relative tolerance for error estimation
     r_tol: f64,
     /// Strategy for updating the Jacobian
@@ -165,7 +166,7 @@ impl<const D: usize> Rosenbrock23<D> {
         Self {
             c32: 6.0 + 2.0_f64.sqrt(),
             d: 1.0 / (2.0 + 2.0_f64.sqrt()),
-            a_tol: 1e-6,
+            a_tol: SVector::<f64, D>::from_element(1e-6),
             r_tol: 1e-3,
             jacobian_strategy: JacobianUpdateStrategy::default(),
             cached_jacobian: None,
@@ -175,18 +176,6 @@ impl<const D: usize> Rosenbrock23<D> {
         }
     }
 
-    /// Set the absolute tolerance
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
-        self.a_tol = a_tol;
-        self
-    }
-
-    /// Set the relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
-        self.r_tol = r_tol;
-        self
-    }
-
     /// Set the Jacobian update strategy
     pub fn jacobian_strategy(mut self, strategy: JacobianUpdateStrategy) -> Self {
         self.jacobian_strategy = strategy;
@@ -209,6 +198,31 @@ impl<const D: usize> Rosenbrock23<D> {
     }
 }
 
+impl<const D: usize> Tolerances<D> for Rosenbrock23<D> {
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
+        self.a_tol = SVector::<f64, D>::from_element(a_tol);
+        self
+    }
+
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+        self.a_tol = a_tol;
+        self
+    }
+
+    fn r_tol(mut self, r_tol: f64) -> Self {
+        self.r_tol = r_tol;
+        self
+    }
+}
+
 impl<const D: usize> Default for Rosenbrock23<D> {
     fn default() -> Self {
         Self::new()
@@ -216,10 +230,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;
@@ -273,8 +298,7 @@ impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
         // Compute scalar error estimate using weighted norm
         let mut error_sum = 0.0;
         for i in 0..D {
-            let scale =
-                self.a_tol + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs());
+            let scale = self.a_tol[i] + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs());
             let weighted_error = error_vec[i].re() / scale;
             error_sum += weighted_error * weighted_error;
         }
@@ -420,12 +444,21 @@ mod tests {
         // So log(error) = 2*log(h) + const
         // Slope should be approximately 2
         for i in 0..errors.len() - 1 {
-            let rate =
-                (errors[i].log10() - errors[i + 1].log10()) / (step_sizes[i].log10() - step_sizes[i + 1].log10());
-            println!("Convergence rate between h={} and h={}: {}", step_sizes[i], step_sizes[i+1], rate);
+            let rate = (errors[i].log10() - errors[i + 1].log10())
+                / (step_sizes[i].log10() - step_sizes[i + 1].log10());
+            println!(
+                "Convergence rate between h={} and h={}: {}",
+                step_sizes[i],
+                step_sizes[i + 1],
+                rate
+            );
 
             // Should be close to 2 for a 2nd order method
-            assert!(rate > 1.8 && rate < 2.2, "Convergence rate {} is not close to 2", rate);
+            assert!(
+                rate > 1.8 && rate < 2.2,
+                "Convergence rate {} is not close to 2",
+                rate
+            );
         }
     }
 
@@ -443,7 +476,7 @@ mod tests {
         let analytical = |t: f64| 0.5 * (1.01 * t).exp();
 
         // Test convergence with Julia's step sizes: (1/2)^(6:-1:3) = [1/64, 1/32, 1/16, 1/8]
-        let step_sizes = vec![1.0/64.0, 1.0/32.0, 1.0/16.0, 1.0/8.0];
+        let step_sizes = vec![1.0 / 64.0, 1.0 / 32.0, 1.0 / 16.0, 1.0 / 8.0];
         let mut errors = Vec::new();
 
         for &h in &step_sizes {
@@ -481,8 +514,11 @@ mod tests {
         println!("Orders per step refinement: {:?}", orders);
 
         // Julia tests: @test sim.𝒪est[:final]≈2 atol=0.2
-        assert!((avg_order - 2.0).abs() < 0.2,
-                "Convergence order {:.2} not within 0.2 of expected order 2", avg_order);
+        assert!(
+            (avg_order - 2.0).abs() < 0.2,
+            "Convergence order {:.2} not within 0.2 of expected order 2",
+            avg_order
+        );
     }
 
     #[test]
@@ -505,11 +541,16 @@ mod tests {
         let mut problem = Problem::new(ode, rb23, controller);
         let solution = problem.solve();
 
-        println!("Adaptive solve completed in {} steps", solution.states.len());
+        println!(
+            "Adaptive solve completed in {} steps",
+            solution.states.len()
+        );
 
-        assert!(solution.states.len() < 20,
-                "Adaptive solve should complete in less than 20 steps, got {}",
-                solution.states.len());
+        assert!(
+            solution.states.len() < 20,
+            "Adaptive solve should complete in less than 20 steps, got {}",
+            solution.states.len()
+        );
 
         // Verify final value is accurate
         let analytical = 0.5 * (1.01_f64 * 1.0).exp();

src/integrator/vern7.rs

@@ -1,9 +1,16 @@
 use nalgebra::SVector;
 
 use super::super::ode::ODE;
-use super::{Integrator, StepResult};
+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)
@@ -73,7 +80,7 @@ pub struct Vern7<const D: usize> {
 impl<const D: usize> Default for Vern7<D>
 where
     Vern7<D>: Integrator<D>,
- {
+{
     fn default() -> Self {
         Self::new()
     }
@@ -92,21 +99,31 @@ where
             r_tol: 1e-8,
         }
     }
+}
 
-    /// Set absolute tolerance (same value for all components)
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for Vern7<D>
+where
+    Vern7<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
         self.a_tol = SVector::<f64, D>::from_element(a_tol);
         self
     }
 
-    /// Set absolute tolerance (different value per component)
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
         self.a_tol = a_tol;
         self
     }
 
-    /// Set relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
+    fn r_tol(mut self, r_tol: f64) -> Self {
         self.r_tol = r_tol;
         self
     }
@@ -124,6 +141,8 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     // Stage 8: [a81, 0, a83, a84, a85, a86, a87]
     // Stage 9: [a91, 0, a93, a94, a95, a96, a97, a98]
     // Stage 10: [a101, 0, a103, a104, a105, a106, a107, 0, 0]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
     const A: &'a [f64] = &[
         // Stage 2
         0.005,
@@ -160,6 +179,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     ];
 
     // Error estimate weights (difference between 7th and 6th order: b - b*)
+    #[rustfmt::skip]
     const B_ERROR: &'a [f64] = &[
         0.002547011879931045,   // b1 - b*1
         0.0,                    // b2 - b*2
@@ -174,6 +194,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     ];
 
     // Time nodes (c coefficients)
+    #[rustfmt::skip]
     const C: &'a [f64] = &[
         0.0,                  // c1
         0.005,                // c2
@@ -193,6 +214,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
 
 impl<'a, const D: usize> Vern7ExtraStages<'a> for Vern7<D> {
     // Time nodes for extra stages
+    #[rustfmt::skip]
     const C_EXTRA: &'a [f64] = &[
         1.0,    // c11
         0.29,   // c12
@@ -279,24 +301,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
@@ -342,6 +380,8 @@ where
         }
     }
 
+    // One interpolation coefficient per line, likewise.
+    #[rustfmt::skip]
     fn interpolate<T: Scalar>(
         &self,
         t_start: T,
@@ -459,7 +499,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];
@@ -553,7 +593,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;
@@ -711,7 +751,7 @@ mod tests {
         }
 
         let y0 = Vector1::new(1.0);
-        let t_end: f64 = 1.0;  // Longer interval to get larger errors
+        let t_end: f64 = 1.0; // Longer interval to get larger errors
         let exact = t_end.exp();
 
         let step_sizes: [f64; 3] = [0.2, 0.1, 0.05];
@@ -738,16 +778,8 @@ mod tests {
         let ratio2 = errors[1] / errors[2];
 
         // Allow some tolerance (expect ratio between 64 and 256)
-        assert!(
-            ratio1 > 64.0 && ratio1 < 256.0,
-            "First ratio: {}",
-            ratio1
-        );
-        assert!(
-            ratio2 > 64.0 && ratio2 < 256.0,
-            "Second ratio: {}",
-            ratio2
-        );
+        assert!(ratio1 > 64.0 && ratio1 < 256.0, "First ratio: {}", ratio1);
+        assert!(ratio2 > 64.0 && ratio2 < 256.0, "Second ratio: {}", ratio2);
     }
 
     #[test]
@@ -795,7 +827,7 @@ mod tests {
             Vector2::new(y[1], -y[0])
         }
 
-        let y0 = Vector2::new(1.0, 0.0);  // Start at maximum displacement, zero velocity
+        let y0 = Vector2::new(1.0, 0.0); // Start at maximum displacement, zero velocity
 
         // Period of harmonic oscillator is 2π
         let period = 2.0 * std::f64::consts::PI;
@@ -821,7 +853,10 @@ mod tests {
 
         println!("Initial energy: {}", energy_initial);
         println!("Final energy: {}", energy_final);
-        println!("Energy drift after {} periods: {:.2e}", num_periods, energy_drift);
+        println!(
+            "Energy drift after {} periods: {:.2e}",
+            num_periods, energy_drift
+        );
         println!("Number of steps: {}", solution.times.len());
 
         // Energy should be conserved to high precision (< 1e-7 relative error over 1000 periods)

src/lib.rs

@@ -16,10 +16,11 @@ 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::integrator::{Integrator, StepResult, Tolerances};
     pub use super::ode::ODE;
-    pub use super::integrator::{Integrator, StepResult};
     pub use super::problem::{Problem, Solution};
     pub use super::scalar::Scalar;
 }
@@ -51,13 +52,13 @@ mod tests {
         let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
         let controller = PIController::default();
 
-        let value_too_high = FnCallback {
+        let mut value_too_high = FnCallback {
             event: &|t: f64, _y: &Vector2<f64>, _p: &Params| 5.0 - t,
             effect: &stop,
         };
 
         // Solve the problem
-        let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high);
+        let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut value_too_high);
         let solution = problem.solve();
 
         // Can interpolate solutions to whatever you want
@@ -321,6 +322,48 @@ mod autodiff_tests {
         }
     }
 
+    /// Derivatives survive a callback that changes the state.
+    ///
+    /// `y' = y` from `y0`, halved the moment it reaches 2, run to `t = 2`. It
+    /// reaches 2 twice in that interval — at `ln 2` and at `2 ln 2` — and `n`
+    /// halvings leave `y(2) = e^2 y0 / 2^n` whatever the times were, so the
+    /// solution and its sensitivity to `y0` are both `e^2 / 4` exactly.
+    ///
+    /// The crossing time is located in real arithmetic and carries no
+    /// derivative. This problem is one where that costs nothing: the state at
+    /// the crossing is fixed at 2 by the event itself, and the final state turns
+    /// out not to depend on when the crossing happened. A problem where the
+    /// event time does matter needs the sensitivity of the event condition,
+    /// which forward mode through Brent's iterations would not give.
+    #[test]
+    fn test_ad_through_a_callback_jump() {
+        fn derivative<T: Scalar>(_t: T, y: Vector1<T>, _p: &()) -> Vector1<T> {
+            y
+        }
+
+        let y0 = Dual64::from_re(1.0).derivative();
+        let mut halve = FnCallback {
+            event: &|_t: Dual64, y: &Vector1<Dual64>, _p: &()| Dual64::lift(2.0) - y[0],
+            effect: &|ode: &mut ODE<'_, Dual64, 1, ()>| ode.y /= Dual64::lift(2.0),
+        };
+
+        let ode = ODE::new(
+            &derivative,
+            Dual64::lift(0.0),
+            Dual64::lift(2.0),
+            Vector1::new(y0),
+            (),
+        );
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let mut problem =
+            Problem::new(ode, vern7, PIController::default()).with_callback(&mut halve);
+        let solution = problem.solve();
+        let y = solution.states.last().unwrap()[0];
+
+        assert_relative_eq!(y.re, 2.0f64.exp() / 4.0, max_relative = 1e-10);
+        assert_relative_eq!(y.eps, 2.0f64.exp() / 4.0, max_relative = 1e-10);
+    }
+
     /// The same, through Vern7 rather than Dormand-Prince, which exercises a
     /// path Dormand-Prince does not have: Vern7 keeps six of its sixteen stages
     /// back and computes them only when something interpolates inside a step.
@@ -333,7 +376,13 @@ mod autodiff_tests {
         }
 
         let y0 = Dual64::from_re(1.0).derivative();
-        let ode = ODE::new(&derivative, Dual64::lift(0.0), Dual64::lift(2.0), Vector1::new(y0), ());
+        let ode = ODE::new(
+            &derivative,
+            Dual64::lift(0.0),
+            Dual64::lift(2.0),
+            Vector1::new(y0),
+            (),
+        );
         let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
         let mut problem = Problem::new(ode, vern7, PIController::default());
         let solution = problem.solve();

src/ode.rs

@@ -23,6 +23,9 @@ pub struct ODE<'a, T: Scalar, const D: usize, P> {
     pub t0: T,
     pub t_end: T,
     pub h: T,
+
+    /// Set when the solution has been ended short of `t_end`, by a callback
+    /// effect calling [`terminate`](ODE::terminate).
     pub finished: bool,
 }
 
@@ -45,6 +48,15 @@ impl<'a, T: Scalar, const D: usize, P> ODE<'a, T, D, P> {
             finished: false,
         }
     }
+
+    /// Ends the solution here, without disturbing `t_end`.
+    ///
+    /// Meant for a callback effect: `Problem::solve` records the state at the
+    /// crossing and then stops. [`stop`](crate::callback::stop) is this as a
+    /// ready-made effect.
+    pub fn terminate(&mut self) {
+        self.finished = true;
+    }
 }
 
 #[cfg(test)]

src/problem.rs

@@ -1,4 +1,11 @@
 //! Assembling an ODE, an integrator and a controller
+//!
+//! Event location is done in real arithmetic: the event condition is a
+//! real-valued question, and Brent's method differentiated through would give
+//! the derivative of its own iteration rather than the sensitivity of the
+//! event, which is wrong in a way that does not announce itself. So the search
+//! runs on `re()` projections and the state at the crossing is produced by one
+//! ordinary step of the integrator, which does carry derivatives.
 
 use nalgebra::SVector;
 use roots::{find_root_brent, SimpleConvergency};
@@ -10,7 +17,22 @@ use super::integrator::Integrator;
 use super::ode::ODE;
 use super::scalar::Scalar;
 
-#[derive(Clone)]
+/// How far into a step the event function of a callback that just fired is read,
+/// as a fraction of the step.
+///
+/// An effect is free to leave its event function sitting on the zero it fired
+/// at — the shadow set swap of a Modified Rodrigues Parameter vector does
+/// exactly that — so reading it at the start of the next step would locate the
+/// same event again, at no advance in time, forever. Reading a little way in
+/// instead gives the sign the event function has after the effect, and the
+/// search for a crossing starts from the same place, so a second crossing
+/// further into the same step is still found.
+///
+/// Small enough that a genuine crossing is very unlikely to hide inside it, and
+/// large enough that the event function has moved clear of its own rounding at
+/// the event.
+const EVENT_PROBE_FRACTION: f64 = 1e-6;
+
 pub struct Problem<'a, T: Scalar, const D: usize, S, P, C>
 where
     S: Integrator<D>,
@@ -19,7 +41,7 @@ where
     ode: ODE<'a, T, D, P>,
     integrator: S,
     controller: C,
-    callbacks: Vec<&'a dyn Callback<T, D, P>>,
+    callbacks: Vec<&'a mut dyn Callback<T, D, P>>,
 }
 
 impl<'a, T: Scalar, const D: usize, S, P, C> Problem<'a, T, D, S, P, C>
@@ -36,6 +58,27 @@ where
         }
     }
 
+    /// Integrates from `t0` to `t_end`, or to wherever a callback ends it.
+    ///
+    /// Each accepted step is checked against every registered callback. When an
+    /// event function changes sign across a step, the step is shortened to land
+    /// on the crossing, the state there is recorded, and the callback's effect
+    /// runs with the ODE at that point — so an effect may change the state, the
+    /// parameters, or end the solution, and the integration carries on from
+    /// whatever it leaves behind.
+    ///
+    /// If more than one event falls inside one step, the earliest fires and the
+    /// step ends there; the others are still ahead and are found again on the
+    /// following step. Ties go to the order the callbacks were registered.
+    ///
+    /// A callback does not fire twice at the same time: an effect may leave its
+    /// event function sitting on the zero it fired at, and the step that begins
+    /// at a firing searches from a hair inside itself rather than from its start.
+    /// A second crossing later in that step is still found — a bouncing ball
+    /// keeps bouncing — but one inside that first millionth of the step is not.
+    ///
+    /// An even number of crossings inside one step is missed, in that the step's
+    /// two ends agree in sign and nothing looks in between.
     pub fn solve(&mut self) -> Solution<'_, T, S, D, P> {
         let mut convergency = SimpleConvergency {
             eps: 1e-12f64,
@@ -44,15 +87,17 @@ where
         let mut times: Vec<T> = vec![self.ode.t];
         let mut states: Vec<SVector<T, D>> = vec![self.ode.y];
         let mut dense_coefficients: Vec<RefCell<Vec<SVector<T, D>>>> = Vec::new();
-        while self.ode.t.re() < self.ode.t_end.re() {
-            if (self.ode.t + self.controller.next_step_guess().extract()).re()
-                > self.ode.t_end.re()
+        // When each callback last fired, which is where its next search for a
+        // crossing starts from rather than at — see `EVENT_PROBE_FRACTION`.
+        let mut last_fired: Vec<Option<f64>> = vec![None; self.callbacks.len()];
+        while !self.ode.finished && self.ode.t.re() < self.ode.t_end.re() {
+            if (self.ode.t + self.controller.next_step_guess().extract()).re() > self.ode.t_end.re()
             {
                 // If the next step would go past the end, then just set it to the end
                 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 (new_y, curr_step, dense_option) = if self.integrator.adaptive() {
                 // First, we try stepping with the "next step guess" to get the error
                 let mut result = self
                     .integrator
@@ -85,32 +130,98 @@ where
                     result.dense,
                 )
             };
-            if !self.callbacks.is_empty() {
-                // Check for events occurring
-                for callback in &self.callbacks {
-                    // The event condition is a real-valued question, so it is
-                    // asked in the reals — see this module's header.
-                    let before = callback.event(self.ode.t, &self.ode.y, &self.ode.params).re();
-                    let after = callback
-                        .event(self.ode.t + curr_step, &new_y, &self.ode.params)
-                        .re();
-                    if before * after < 0.0 {
-                        // If the event crossed zero, then find the root
-                        let f = |test_h: f64| {
-                            let h = T::lift(test_h);
-                            let test_y = self.integrator.step(&self.ode, h).y;
-                            callback.event(self.ode.t + h, &test_y, &self.ode.params).re()
-                        };
-                        let root =
-                            find_root_brent(0.0f64, curr_step.re(), &f, &mut convergency).unwrap();
-                        curr_step = T::lift(root);
-                        let result = self.integrator.step(&self.ode, curr_step);
-                        new_y = result.y;
-                        dense_option = result.dense;
-                        callback.effect(&mut self.ode);
+            // Look for the earliest event in the step. Nothing is committed
+            // while looking, because an event function is asked about the state
+            // at the start of the step as well as at the end.
+            let mut earliest: Option<(usize, f64)> = None;
+            for (index, (callback, fired)) in
+                self.callbacks.iter().zip(last_fired.iter()).enumerate()
+            {
+                // A callback that fired at this step's start searches from a
+                // little way in instead, so that it cannot find the event it
+                // just handled again.
+                let search_start = if *fired == Some(self.ode.t.re()) {
+                    curr_step * EVENT_PROBE_FRACTION
+                } else {
+                    T::lift(0.0)
+                };
+                // The event condition is a real-valued question, so it is
+                // asked in the reals — see this module's header.
+                let before = if search_start.re() == 0.0 {
+                    callback
+                        .event(self.ode.t, &self.ode.y, &self.ode.params)
+                        .re()
+                } else {
+                    let probe = self.integrator.step(&self.ode, search_start).y;
+                    callback
+                        .event(self.ode.t + search_start, &probe, &self.ode.params)
+                        .re()
+                };
+                let after = callback
+                    .event(self.ode.t + curr_step, &new_y, &self.ode.params)
+                    .re();
+                let crossing = if before == 0.0 {
+                    // Sitting exactly on the event. At the start of a step that
+                    // is a crossing, and it fires there; a little way into one
+                    // it is an event function this step can do nothing with,
+                    // because the effect has already run at its zero.
+                    (search_start.re() == 0.0).then_some(0.0)
+                } else if after == 0.0 {
+                    Some(curr_step.re())
+                } else if before * after < 0.0 {
+                    let f = |test_h: f64| {
+                        let h = T::lift(test_h);
+                        let test_y = self.integrator.step(&self.ode, h).y;
+                        callback
+                            .event(self.ode.t + h, &test_y, &self.ode.params)
+                            .re()
+                    };
+                    Some(
+                        find_root_brent(search_start.re(), curr_step.re(), &f, &mut convergency)
+                            .unwrap(),
+                    )
+                } else {
+                    None
+                };
+                // Two events in one step: the earlier one fires and the step
+                // ends there, so the later one is still ahead and is found
+                // again on the next step. A tie goes to registration order.
+                if let Some(h_event) = crossing {
+                    if earliest.is_none_or(|(_, h)| h_event < h) {
+                        earliest = Some((index, h_event));
                     }
                 }
             }
+
+            if let Some((index, h_event)) = earliest {
+                // Land on the crossing and record it, before the effect runs:
+                // the trajectory does reach this state, and a terminal effect
+                // should leave it as the last point of the solution.
+                if h_event != 0.0 {
+                    let result = self.integrator.step(&self.ode, T::lift(h_event));
+                    self.ode.t += T::lift(h_event);
+                    self.ode.y = result.y;
+                    times.push(self.ode.t);
+                    states.push(self.ode.y);
+                    dense_coefficients.push(RefCell::new(result.dense.unwrap()));
+                }
+                last_fired[index] = Some(self.ode.t.re());
+                let before_effect = self.ode.y;
+                self.callbacks[index].effect(&mut self.ode);
+                if self.ode.y != before_effect {
+                    // The state jumped. Both sides of it are part of the
+                    // solution, so both are recorded at the same time, with a
+                    // zero-length interval between them. The controller's error
+                    // history describes the solution up to the jump and says
+                    // nothing about the one after it.
+                    times.push(self.ode.t);
+                    states.push(self.ode.y);
+                    dense_coefficients.push(RefCell::new(vec![before_effect, self.ode.y]));
+                    self.controller.reset();
+                }
+                continue;
+            }
+
             self.ode.y = new_y;
             self.ode.t += curr_step;
             times.push(self.ode.t);
@@ -127,12 +238,30 @@ where
         }
     }
 
-    pub fn with_callback(mut self, callback: &'a dyn Callback<T, D, P>) -> Self {
+    /// Registers a callback.
+    ///
+    /// Callbacks are watched in the order they are registered, which decides
+    /// which one fires when two events land at the same time.
+    pub fn with_callback(mut self, callback: &'a mut dyn Callback<T, D, P>) -> Self {
         self.callbacks.push(callback);
         self
     }
 }
 
+/// What an integration produced: the steps it took, and the means to evaluate
+/// between them.
+///
+/// `times` and `states` run in step order and are the same length. `dense` has
+/// one fewer entry, the coefficients for interpolating inside the step from
+/// `times[i]` to `times[i + 1]`.
+///
+/// A callback effect that changes the state puts the same time in `times`
+/// twice, with the state before the jump and the state after it. The pair is
+/// what the trajectory does: it reaches one value and leaves from the other.
+/// [`interpolate`](Solution::interpolate) at exactly that time reports the
+/// second of them, the state the effect left, since that is the one the
+/// trajectory carries forward; the value it arrived with is the entry before it
+/// in `states`.
 pub struct Solution<'a, T: Scalar, S, const D: usize, P>
 where
     S: Integrator<D>,
@@ -148,6 +277,10 @@ impl<T: Scalar, S, const D: usize, P> Solution<'_, T, S, D, P>
 where
     S: Integrator<D>,
 {
+    /// The state at `t`, interpolated inside whichever step contains it.
+    ///
+    /// At the time of a jump this is the state the callback's effect left. Panics
+    /// if `t` is outside the solution.
     pub fn interpolate(&self, t: T) -> SVector<T, D> {
         // First check that the t is within bounds
         let last = self.times.last().unwrap().re();
@@ -164,7 +297,16 @@ where
 
         // Then find the two t values closest to the desired t
         match times.binary_search_by(|x| x.total_cmp(&t.re())) {
-            Ok(index) => self.states[index],
+            Ok(index) => {
+                // A jump records its time twice and the search may land on
+                // either entry. Walk to the last of them, which is the state the
+                // effect left and the one the trajectory goes on with.
+                let mut index = index;
+                while index + 1 < times.len() && times[index + 1] == times[index] {
+                    index += 1;
+                }
+                self.states[index]
+            }
             Err(end_index) => {
                 // Then send that to the integrator
                 let t_start = self.times[end_index - 1];
@@ -172,16 +314,24 @@ where
                 let y_start = self.states[end_index - 1];
                 let h = t_end - t_start;
 
+                if h.re() == 0.0 {
+                    // The two sides of a jump. Any `t` strictly inside this
+                    // interval would have to be the jump's own time, which the
+                    // search above answers exactly, so this is only reached if
+                    // a step of zero length was recorded some other way.
+                    return y_start;
+                }
+
                 // 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 +339,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(
@@ -219,8 +369,19 @@ mod tests {
     use crate::callback::{stop, FnCallback};
     use crate::controller::PIController;
     use crate::integrator::dormand_prince::DormandPrince45;
+    use crate::integrator::vern7::Vern7;
+    use crate::integrator::Tolerances;
     use approx::assert_relative_eq;
-    use nalgebra::Vector3;
+    use nalgebra::{Matrix3, Vector1, Vector2, Vector3};
+    use std::cell::Cell;
+    use std::f64::consts::PI;
+    use std::rc::Rc;
+
+    /// `dy/dt = 1`, so the state is the time elapsed and an event on the state
+    /// is an event on the clock.
+    fn unit_rate(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+        Vector1::new(1.0)
+    }
 
     #[test]
     fn test_problem() {
@@ -258,12 +419,12 @@ mod tests {
         let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-5);
         let controller = PIController::default();
 
-        let value_too_high = FnCallback {
+        let mut value_too_high = FnCallback {
             event: &|_: f64, y: &SVector<f64, 3>, _: &Params| 10.0 - y[0],
             effect: &stop,
         };
 
-        let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high);
+        let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut value_too_high);
         let solution = problem.solve();
 
         assert_relative_eq!(
@@ -273,6 +434,415 @@ mod tests {
         );
     }
 
+    /// An effect that writes to `ode.y` is what the integration carries on from.
+    #[test]
+    fn test_effect_changes_the_state() {
+        let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10);
+        let mut jump = FnCallback {
+            event: &|_t: f64, y: &Vector1<f64>, _p: &()| 5.0 - y[0],
+            effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
+        };
+
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
+        let solution = problem.solve();
+
+        // Five units of integration before the jump, five after.
+        assert_relative_eq!(solution.times.last().unwrap(), &10.0, max_relative = 1e-12);
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            105.0,
+            max_relative = 1e-9
+        );
+    }
+
+    /// A jump is recorded as two states at the same time, in the order the
+    /// trajectory visits them.
+    #[test]
+    fn test_jump_records_both_sides() {
+        let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10);
+        let mut jump = FnCallback {
+            event: &|_t: f64, y: &Vector1<f64>, _p: &()| 5.0 - y[0],
+            effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
+        };
+
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
+        let solution = problem.solve();
+
+        let at_crossing: Vec<usize> = (0..solution.times.len())
+            .filter(|&i| (solution.times[i] - 5.0).abs() < 1e-9)
+            .collect();
+        assert_eq!(at_crossing.len(), 2, "the two sides of the jump");
+        let (left, right) = (at_crossing[0], at_crossing[1]);
+        assert_eq!(right, left + 1);
+        assert_relative_eq!(solution.states[left][0], 5.0, max_relative = 1e-9);
+        assert_relative_eq!(solution.states[right][0], 100.0, max_relative = 1e-12);
+
+        // Interpolation either side of the jump follows the branch it is on,
+        // and at the jump itself reports the state the effect left.
+        assert_relative_eq!(solution.interpolate(4.9)[0], 4.9, max_relative = 1e-9);
+        assert_relative_eq!(solution.interpolate(5.1)[0], 100.1, max_relative = 1e-9);
+        assert_relative_eq!(
+            solution.interpolate(solution.times[left])[0],
+            100.0,
+            max_relative = 1e-12
+        );
+    }
+
+    /// `stop` ends the solution at the crossing, and leaves `t_end` alone.
+    #[test]
+    fn test_stop_ends_at_the_crossing() {
+        let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10);
+        let mut stopper = FnCallback {
+            event: &|_t: f64, y: &Vector1<f64>, _p: &()| 5.0 - y[0],
+            effect: &stop,
+        };
+
+        let mut problem =
+            Problem::new(ode, dp45, PIController::default()).with_callback(&mut stopper);
+        let solution = problem.solve();
+
+        assert_relative_eq!(solution.times.last().unwrap(), &5.0, max_relative = 1e-12);
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            5.0,
+            max_relative = 1e-12
+        );
+        // Nothing past the event, and the problem's own end time is unchanged.
+        assert!(solution.times.iter().all(|&t| t <= 5.0 + 1e-12));
+        assert_eq!(solution.ode.t_end, 10.0);
+        assert!(solution.ode.finished);
+    }
+
+    /// `effect` takes `&mut self`, so a callback can count its own firings and
+    /// act only on the first one.
+    #[test]
+    fn test_callback_can_remember_firing() {
+        struct FireOnce {
+            fired: usize,
+        }
+
+        impl Callback<f64, 1, ()> for FireOnce {
+            fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
+                // Crosses at t = 1, 3, 5 and 7.
+                (t * PI / 2.0).cos()
+            }
+
+            fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) {
+                self.fired += 1;
+                if self.fired == 1 {
+                    ode.y = Vector1::new(100.0);
+                }
+            }
+        }
+
+        let ode = ODE::new(&unit_rate, 0.0, 9.0, Vector1::new(0.0), ());
+        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10);
+        let mut once = FireOnce { fired: 0 };
+        // `dy/dt = 1` has no local error at all, so the step size runs away to
+        // `h_max`; capped at half a period, since a step long enough to contain
+        // two crossings sees no sign change and finds neither.
+        let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 1.0, 0.9, 1e-4);
+
+        let mut problem = Problem::new(ode, dp45, controller).with_callback(&mut once);
+        let solution = problem.solve();
+        let last = solution.states.last().unwrap()[0];
+
+        assert_eq!(once.fired, 4);
+        // Only the first firing moved the state: 1 unit before the jump, 8 after.
+        assert_relative_eq!(last, 108.0, max_relative = 1e-9);
+    }
+
+    /// An event function that is already zero at the start of a step crosses
+    /// there, and fires there, rather than being passed over.
+    #[test]
+    fn test_event_already_zero_fires_at_the_step_start() {
+        // y starts at zero and the event is y, so the event is met at t = 0.
+        let ode = ODE::new(&unit_rate, 0.0, 1.0, Vector1::new(0.0), ());
+        let mut at_zero = FnCallback {
+            event: &|_t: f64, y: &Vector1<f64>, _p: &()| y[0],
+            effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
+        };
+
+        let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+            .with_callback(&mut at_zero);
+        let solution = problem.solve();
+
+        assert_eq!(solution.times[0], 0.0);
+        assert_eq!(solution.states[0][0], 0.0);
+        assert_eq!(solution.times[1], 0.0);
+        assert_eq!(solution.states[1][0], 100.0);
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            101.0,
+            max_relative = 1e-9
+        );
+    }
+
+    /// Two events inside one step: the earlier one fires, and the later one is
+    /// found again on the step after it. Which callback was registered first
+    /// does not change the order they fire in.
+    #[test]
+    fn test_earliest_event_in_a_step_fires_first() {
+        fn order_of_firings(three_first: bool) -> Vec<f64> {
+            let fired = Rc::new(Cell::new([0.0f64; 2]));
+
+            let at_three = {
+                let fired = Rc::clone(&fired);
+                move |ode: &mut ODE<'_, f64, 1, ()>| {
+                    let mut seen = fired.get();
+                    seen[0] = ode.t;
+                    fired.set(seen);
+                }
+            };
+            let at_five = {
+                let fired = Rc::clone(&fired);
+                move |ode: &mut ODE<'_, f64, 1, ()>| {
+                    let mut seen = fired.get();
+                    seen[1] = ode.t;
+                    fired.set(seen);
+                }
+            };
+
+            let mut three = FnCallback {
+                event: &|t: f64, _y: &Vector1<f64>, _p: &()| 3.0 - t,
+                effect: &at_three,
+            };
+            let mut five = FnCallback {
+                event: &|t: f64, _y: &Vector1<f64>, _p: &()| 5.0 - t,
+                effect: &at_five,
+            };
+
+            let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+            // One step of 10 would swallow both events whole.
+            let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 100000.0, 0.9, 10.0);
+            let mut problem = if three_first {
+                Problem::new(ode, DormandPrince45::new(), controller)
+                    .with_callback(&mut three)
+                    .with_callback(&mut five)
+            } else {
+                Problem::new(ode, DormandPrince45::new(), controller)
+                    .with_callback(&mut five)
+                    .with_callback(&mut three)
+            };
+            problem.solve();
+            fired.get().to_vec()
+        }
+
+        for three_first in [true, false] {
+            let fired = order_of_firings(three_first);
+            assert_relative_eq!(fired[0], 3.0, max_relative = 1e-9);
+            assert_relative_eq!(fired[1], 5.0, max_relative = 1e-9);
+        }
+    }
+
+    /// An effect that leaves the event function sitting on zero does not
+    /// re-trigger the callback at the same point.
+    #[test]
+    fn test_effect_landing_on_its_own_event_terminates() {
+        // The event is |y| - 1 and the effect negates y, so the event value is
+        // unchanged by the effect: the naive reading is that it has crossed
+        // again and should fire again, at no advance in time.
+        let ode = ODE::new(&unit_rate, 0.0, 3.0, Vector1::new(0.0), ());
+        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-10);
+        let mut flip = FnCallback {
+            event: &|_t: f64, y: &Vector1<f64>, _p: &()| y[0].abs() - 1.0,
+            effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = -ode.y,
+        };
+
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
+        let solution = problem.solve();
+
+        // y rises to 1 at t = 1 and flips to -1, which the effect leaves on the
+        // event; it rises back to 1 at t = 3 and flips again. Two jumps, neither
+        // of them repeated at the time it happened.
+        let jumps: Vec<f64> = (1..solution.times.len())
+            .filter(|&i| solution.times[i] == solution.times[i - 1])
+            .map(|i| solution.times[i])
+            .collect();
+        assert_eq!(jumps.len(), 2);
+        assert_relative_eq!(jumps[0], 1.0, max_relative = 1e-9);
+        assert_relative_eq!(jumps[1], 3.0, max_relative = 1e-9);
+        assert_relative_eq!(solution.times.last().unwrap(), &3.0, max_relative = 1e-12);
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            -1.0,
+            max_relative = 1e-8
+        );
+    }
+
+    /// The controller's error history is discarded when an effect changes the
+    /// state, and left alone when it does not.
+    #[test]
+    fn test_controller_is_reset_across_a_jump() {
+        /// Counts the resets `solve` asks for; a `PIController` otherwise.
+        struct Spy {
+            inner: PIController<f64>,
+            resets: Rc<Cell<usize>>,
+        }
+
+        impl Controller<f64> for Spy {
+            fn determine_step(&mut self, h: f64, err: f64) -> TryStep<f64> {
+                self.inner.determine_step(h, err)
+            }
+
+            fn next_step_guess(&self) -> TryStep<f64> {
+                Controller::next_step_guess(&self.inner)
+            }
+
+            fn set_next_step_guess(&mut self, guess: TryStep<f64>) {
+                self.inner.set_next_step_guess(guess)
+            }
+
+            fn reset(&mut self) {
+                self.resets.set(self.resets.get() + 1);
+                Controller::reset(&mut self.inner);
+            }
+        }
+
+        fn resets_during(effect: &dyn Fn(&mut ODE<'_, f64, 1, ()>)) -> usize {
+            let resets = Rc::new(Cell::new(0));
+            let spy = Spy {
+                inner: PIController::default(),
+                resets: Rc::clone(&resets),
+            };
+            let mut callback = FnCallback {
+                event: &|_t: f64, y: &Vector1<f64>, _p: &()| 5.0 - y[0],
+                effect,
+            };
+            let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+            let mut problem =
+                Problem::new(ode, DormandPrince45::new(), spy).with_callback(&mut callback);
+            problem.solve();
+            resets.get()
+        }
+
+        assert_eq!(
+            resets_during(&|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0)),
+            1
+        );
+        // An effect that changes nothing leaves a smooth solution behind it.
+        assert_eq!(resets_during(&|_ode: &mut ODE<'_, f64, 1, ()>| {}), 0);
+    }
+
+    /// A bouncing ball, which is the standard exercise of repeated jumps: an
+    /// event every time the ball reaches the floor, and an effect that reverses
+    /// and shrinks its velocity.
+    ///
+    /// The bounce times are known. A ball dropped from `h0` first lands at
+    /// `sqrt(2 h0 / g)` with speed `v = g t`, and each bounce after that is a
+    /// flight of `2 v / g` at the restitution-scaled speed.
+    #[test]
+    fn test_bouncing_ball() {
+        const G: f64 = 9.81;
+        const RESTITUTION: f64 = 0.8;
+
+        fn free_fall(_t: f64, y: Vector2<f64>, g: &f64) -> Vector2<f64> {
+            Vector2::new(y[1], -g)
+        }
+
+        let ode = ODE::new(&free_fall, 0.0, 2.0, Vector2::new(1.0, 0.0), G);
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let mut bounce = FnCallback {
+            event: &|_t: f64, y: &Vector2<f64>, _p: &f64| y[0],
+            effect: &|ode: &mut ODE<'_, f64, 2, f64>| {
+                ode.y = Vector2::new(0.0, -RESTITUTION * ode.y[1]);
+            },
+        };
+
+        let mut problem =
+            Problem::new(ode, vern7, PIController::default()).with_callback(&mut bounce);
+        let solution = problem.solve();
+
+        // Each bounce is a jump, so it is a repeated time in the solution.
+        let bounces: Vec<(f64, f64)> = (1..solution.times.len())
+            .filter(|&i| solution.times[i] == solution.times[i - 1])
+            .map(|i| (solution.times[i], solution.states[i][1]))
+            .collect();
+
+        // Three of them land before t = 2, and the fourth does not.
+        let mut expected_t = (2.0 / G).sqrt();
+        let mut expected_v = RESTITUTION * G * expected_t;
+        assert_eq!(bounces.len(), 3);
+        for (t, v) in &bounces {
+            assert_relative_eq!(*t, expected_t, max_relative = 1e-9);
+            assert_relative_eq!(*v, expected_v, max_relative = 1e-9);
+            expected_t += 2.0 * expected_v / G;
+            expected_v *= RESTITUTION;
+        }
+
+        // Interpolating inside an arc gives free fall from the bounce that
+        // started it, so the dense output survived the jumps.
+        for &(t, v) in &bounces[..2] {
+            for dt in [1e-3, 0.05, 0.1] {
+                assert_relative_eq!(
+                    solution.interpolate(t + dt)[0],
+                    v * dt - 0.5 * G * dt * dt,
+                    epsilon = 1e-9
+                );
+            }
+        }
+    }
+
+    /// The case the whole feature is for: a rotating rigid body parameterized by
+    /// Modified Rodrigues Parameters, whose parameters are swapped for the
+    /// shadow set before they run into the singularity at a full turn.
+    ///
+    /// For a constant angular velocity about one axis the answer is known.
+    /// Starting from zero, the MRP vector is `tan(theta / 4)` along the rotation
+    /// axis, with `theta = |omega| t`, so it reaches magnitude one at half a
+    /// turn. The shadow set `-sigma / |sigma|^2` describes the same attitude and
+    /// obeys the same kinematics, and past half a turn it is the branch that
+    /// stays bounded: `-cot(theta / 4)` along the axis.
+    #[test]
+    fn test_mrp_shadow_set_swap() {
+        /// MRP kinematics, `sigma' = B(sigma) omega / 4`.
+        fn mrp_rate(_t: f64, s: Vector3<f64>, omega: &Vector3<f64>) -> Vector3<f64> {
+            let skew = Matrix3::new(0.0, -s[2], s[1], s[2], 0.0, -s[0], -s[1], s[0], 0.0);
+            let b = Matrix3::identity() * (1.0 - s.dot(&s)) + skew * 2.0 + s * s.transpose() * 2.0;
+            b * omega * 0.25
+        }
+
+        let omega = Vector3::new(0.0, 0.0, 1.0);
+        let t_end = 5.0; // past half a turn, short of a full one
+        let ode = ODE::new(&mrp_rate, 0.0, t_end, Vector3::zeros(), omega);
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+
+        let mut shadow = FnCallback {
+            event: &|_t: f64, s: &Vector3<f64>, _p: &Vector3<f64>| 1.0 - s.norm(),
+            effect: &|ode: &mut ODE<'_, f64, 3, Vector3<f64>>| {
+                ode.y = -ode.y / ode.y.norm_squared();
+            },
+        };
+        let mut problem =
+            Problem::new(ode, vern7, PIController::default()).with_callback(&mut shadow);
+        let solution = problem.solve();
+
+        // The swap happened at half a turn, and only once. A jump is the one
+        // thing that puts the same time in the solution twice.
+        let jumps: Vec<f64> = (1..solution.times.len())
+            .filter(|&i| solution.times[i] == solution.times[i - 1])
+            .map(|i| solution.times[i])
+            .collect();
+        assert_eq!(jumps.len(), 1);
+        assert_relative_eq!(jumps[0], PI, max_relative = 1e-9);
+
+        // The parameters stayed away from the singularity throughout.
+        assert!(solution.states.iter().all(|s| s.norm() <= 1.0 + 1e-9));
+
+        // And the attitude at the end is the one the body is actually in.
+        let expected = -Vector3::new(0.0, 0.0, 1.0 / (t_end / 4.0).tan());
+        for i in 0..3 {
+            assert_relative_eq!(
+                solution.states.last().unwrap()[i],
+                expected[i],
+                epsilon = 1e-9
+            );
+        }
+    }
+
     #[test]
     fn test_with_interpolation() {
         type Params = ();

src/scalar.rs

@@ -115,5 +115,4 @@ mod tests {
         // d/dx x^3 = 3x^2 = 36.75 at x = 3.5
         assert_eq!(x.powi(3).eps, 3.0 * 3.5 * 3.5);
     }
-
 }