Merge branch 'release/0.8.0'

Scalar now inherits num-dual's DualNum<Primitive = f64> and
DualStruct<Real = f64> and declares only lift. re comes from DualStruct.
A blanket impl replaces impl Scalar for f64 and the impl_scalar! macro
over the dual types, so Scalar covers f64 and every num-dual dual type
and a sibling crate built on num-dual can share one bound.

num-dual is a plain dependency rather than optional.

Breaks: the autodiff feature is removed, so a dependent naming it in its
feature list no longer resolves. Scalar's blanket impl makes any manual
impl of it a conflicting implementation, and Scalar no longer declares
re, so a qualified Scalar::re(&x) call must become x.re().

Cargo.toml

@@ -1,29 +1,25 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.3.1"
+version = "0.8.0"
 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"]
-autodiff = ["dep:num-dual"]
 
 [dependencies]
 serde = { version = "1.0", features = ["derive"] }
 nalgebra = { version = "0.35", features = ["serde-serialize"] }
 num-traits = "0.2.19"
-num-dual = { version = "0.15", optional = true }
+num-dual = "0.15"
 roots = "0.0.8"
 
 [dev-dependencies]
 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,25 +72,159 @@ 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.
+
+### Borrowed or owned
+
+`with_callback` borrows, which is what lets you read a stateful callback's own fields once the solve
+is done:
+
+```rust
+# use nalgebra::Vector1;
+# use ordinary_diffeq::prelude::*;
+struct Count {
+    fired: usize,
+}
+
+impl Callback<f64, 1, ()> for Count {
+    fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
+        1.0 - t
+    }
+
+    fn effect(&mut self, _ode: &mut ODE<'_, f64, 1, ()>) {
+        self.fired += 1;
+    }
+}
+
+fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+    Vector1::new(1.0)
+}
+
+let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(0.0), ());
+let mut counter = Count { fired: 0 };
+let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+    .with_callback(&mut counter);
+problem.solve();
+
+assert_eq!(counter.fired, 1);
+```
+
+`with_owned_callback` takes a `Box` instead, for a caller that has nothing to lend. A type that wraps
+this one has to keep the `Problem` in a field, since `solve` hands back a solution borrowing from it,
+and it then cannot also hold a callback in a field and lend it to the problem — that is a
+self-referential struct, and moving the wrapper would dangle the reference. Moving the callback into
+the problem sidesteps it:
+
+```rust,ignore
+let problem = Problem::new(ode, integrator, controller)
+    .with_owned_callback(Box::new(my_callback));
+self.problem = Some(problem);
+self.problem.as_mut().unwrap().solve()
+```
+
+An owned callback may hold borrows of its own, so a boxed `FnCallback` works. Nothing can read its
+state afterwards; a callback with something to report either stays borrowed, or shares an
+`Rc<Cell<_>>` with its owner, which can be read while the solution is still alive.
+
+## 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,
+State, time and step size are generic over the `Scalar` trait, which is a blanket extension of
-and so do the forward-mode dual numbers of [`num-dual`](https://crates.io/crates/num-dual) (the
+`DualNum` and `DualStruct` from [`num-dual`](https://crates.io/crates/num-dual): those two supply
-`autodiff` feature, on by default). Solving at a dual scalar propagates exact derivatives alongside
+the arithmetic, the transcendentals and `re`, which discards derivative information, and `Scalar`
-the state.
+adds only `lift`. Both are implemented for `f64` as well as for every one of that crate's
+forward-mode dual numbers, so `Scalar` covers all of them, and a crate depending on this one
+alongside another built on `num-dual` can satisfy both with one bound. Solving at a dual scalar
+propagates exact derivatives alongside the state.
 
 The only change needed is to the signature of the right-hand side. Here is the pendulum from above,
 before and after:

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,41 @@ 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
+  - extended in 0.7.0 so a `Problem` can own a callback as well as borrow one,
+    which is what a type that stores the problem needs
+  - 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 +108,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 +126,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 +138,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 +153,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 +161,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 +169,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 +209,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 +227,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 +273,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 +293,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 +370,20 @@ Each algorithm implementation should include:
 
 ## Progress Tracking
 
-Total Features: 38
+Total Features: 42
-- Tier 1: 8 features (4/8 complete) ✅
+- Tier 1: 11 features (7/11 complete) ✅
-- Tier 2: 12 features (0/12 complete)
+- Tier 2: 11 features (0/11 complete)
-- Tier 3: 18 features (0/18 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,333 @@
+# Feature: Callback Effects That Can Change the State
+
+**Status**: completed (2026-10-01), extended with owned callbacks (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.
+
+## Follow-up: a problem can own a callback (0.7.0)
+
+`Problem` held `Vec<&'a mut dyn Callback>`, so something else had to own each
+callback and outlive the problem. That is fine from a function body and impossible
+from a type that stores the problem — and a wrapper must store it, because
+`solve(&mut self) -> Solution<'_, ..>` hands back a solution that borrows from it,
+so the problem has to live somewhere that outlives the solution, and a field is
+the only such place. Such a wrapper could not then register a callback of its own:
+
+```rust
+struct Wrapper<'a, ..> {
+    my_callback: MyCallback,           // owned here
+    problem: Option<Problem<'a, ..>>,  // and borrowed from, here
+}
+```
+
+That is self-referential, and moving the wrapper would dangle the reference, so
+the compiler is right to reject it. The workaround in the dependent crate was a
+zero-sized callback, leaked to get a `&'static mut`.
+
+`with_owned_callback(Box<dyn Callback<T, D, P> + 'o>)` is the addition. Storage is
+one list of either kind:
+
+```rust
+enum Registered<'a, 'o, T: Scalar, const D: usize, P> {
+    Borrowed(&'a mut dyn Callback<T, D, P>),
+    Owned(Box<dyn Callback<T, D, P> + 'o>),
+}
+```
+
+`Registered` implements `Callback` itself, forwarding both methods, which is why
+none of the event handling in `solve` changed: the immutable scan for the earliest
+crossing goes through `event` and the single effect that fires goes through
+`effect`, whichever kind is holding it. Root finding, jump recording, the
+controller reset and `EVENT_PROBE_FRACTION` are untouched.
+
+### The second lifetime, which was not in the plan
+
+`Problem` gained a lifetime parameter: `Problem<'a, 'o, T, D, S, P, C>`. This was
+meant to be purely additive and it cannot be, for a reason worth writing down.
+
+Dropping a `Box<dyn Callback>` runs a destructor the compiler cannot see into, so
+dropck requires everything that box's lifetime covers to be alive when the problem
+is dropped. With one lifetime for both kinds, that requirement lands on the
+borrowed callbacks too — and it is type-based, so it applies even to a problem that
+owns no callback at all. The first attempt, a single `'a`, failed exactly here:
+
+```
+error[E0502]: cannot borrow `once.fired` as immutable because it is also borrowed as mutable
+  |
+  | mutable borrow might be used here, when `problem` is dropped
+```
+
+That is `test_callback_can_remember_firing`, reading a borrowed callback's own
+state after the solve, which is the only way a caller learns what a stateful
+callback did and the whole reason the borrowed form exists. Keeping it means the
+owned box needs a lifetime of its own, so that only `'o` has to outlive the drop
+and NLL can still end a borrowed callback's borrow at the problem's last use.
+
+The cost is arity: code that names the type explicitly needs
+`Problem<'a, '_, ..>`. Both lifetimes are inferred, so nothing that only calls
+`Problem::new(..).with_callback(..).solve()` is affected.
+
+### Not done: getting owned callbacks back
+
+Deliberately skipped. Downcasting a `Box<dyn Callback>` to its concrete type needs
+`Any`, which means a `'static` bound and an `as_any` method on every implementor —
+a cost paid by all callbacks for a need that has two cheaper answers already: stay
+borrowed, or share an `Rc<Cell<_>>`, which can be read while the solution is still
+alive.
+
+### Measurement
+
+The benches register no callbacks, so what they exercise is the list's element type
+changing from a reference to a two-variant enum, on a list that is empty. Measured
+by alternating old and new rather than running one after the other, because the
+first attempt at it reported a 4% improvement, which nothing in this change could
+produce: between-sitting drift on this machine is larger than the effect being
+looked for.
+
+| bench | old | new |
+|---|---|---|
+| `bench_orbit` | 356.06, 358.91 us | 357.98, 355.18 us |
+| `bench_simple_1d` | 877.7, 881.1 ns | 877.4, 895.7 ns |
+| `bench_interpolation_1d` | 2.2683, 2.2104 us | 2.2444, 2.2496 us |
+
+Two runs of the same binary differ by up to 2.6% here, and the old and new ranges
+overlap on all three. No measurable change.
+
+### Testing
+
+- an owned callback producing the same times and states as the identical borrowed
+  one, element for element;
+- an owned stateful callback escalating on its second firing, read off the
+  trajectory since its state cannot be read directly;
+- owned and borrowed callbacks on one problem with their events at the same time,
+  where the final state distinguishes registration order, and the borrowed one's
+  count is read after the solve;
+- a boxed `FnCallback`, which holds two closure references, confirming the `+ 'o`
+  rather than `'static`;
+- the wrapper from the top of this section, compiling and solving.
+
+## 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
+/// step and looks for a sign change. When it finds one it shortens the step to
-/// and applies [`effect`](Callback::effect) to the ODE there.
+/// 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(
+            TryStep::NotYetAccepted(prev_step / self.factor_c1.min(factor_11 / self.safety_factor))
-                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 p_term: f64, // Proportional: reaction to current error
-    pub i_term: f64,  // Integral: reaction to error history
+    pub i_term: f64, // Integral: reaction to error history
-    pub d_term: f64,  // Derivative: reaction to error rate of change
+    pub d_term: f64, // Derivative: reaction to error rate of change
 
     // Constraints
-    pub factor_c1: f64,      // 1/min_factor (maximum step decrease)
+    pub factor_c1: f64,     // 1/min_factor (maximum step decrease)
-    pub factor_c2: f64,      // 1/max_factor (maximum step increase)
+    pub factor_c2: f64,     // 1/max_factor (maximum step increase)
-    pub h_max: f64,          // Maximum allowed step size
+    pub h_max: f64,         // Maximum allowed step size
-    pub safety_factor: f64,  // Safety factor (typically 0.9)
+    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_old: f64,   // ε_{n-1}: previous step error
-    pub err_older: f64,   // ε_{n-2}: error two steps ago
+    pub err_older: f64, // ε_{n-2}: error two steps ago
-    pub h_old: f64,       // h_{n-1}: previous step size
+    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_c2
-                self.factor_c1.min(factor_11 / self.safety_factor)
+                .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_old: 0.0,   // No history initially
-            err_older: 0.0,    // No history initially
+            err_older: 0.0, // No history initially
-            h_old: 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.49 / k_plus_1, // p_term: proportional
-            0.34 / k_plus_1,  // i_term: integral
+            0.34 / k_plus_1, // i_term: integral
-            0.10 / k_plus_1,  // d_term: derivative
+            0.10 / k_plus_1, // d_term: derivative
-            10.0,              // max_factor
+            10.0,            // max_factor
-            0.2,               // min_factor
+            0.2,             // min_factor
-            100000.0,          // h_max
+            100000.0,        // h_max
-            0.9,               // safety_factor
+            0.9,             // safety_factor
-            1e-4,              // initial_h
+            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.07,     // p_term (proportional)
-            0.04,      // i_term (integral)
+            0.04,     // i_term (integral)
-            0.01,      // d_term (derivative)
+            0.01,     // d_term (derivative)
-            10.0,      // max_factor
+            10.0,     // max_factor
-            0.2,       // min_factor
+            0.2,      // min_factor
-            100000.0,  // h_max
+            100000.0, // h_max
-            0.9,       // safety_factor
+            0.9,      // safety_factor
-            1e-4,      // initial_h
+            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
+        let pi_variance: f64 =
-            .iter()
+            pi_ratios.iter().map(|r| (r - pi_mean).powi(2)).sum::<f64>() / pi_ratios.len() as f64;
-            .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)
+impl<const D: usize> Tolerances<D> for BS3<D>
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+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)
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
         self.a_tol = a_tol;
         self
     }
 
-    /// Set relative tolerance
+    fn r_tol(mut self, r_tol: f64) -> Self {
-    pub 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;
+    fn order(&self) -> usize {
-    const STAGES: usize = 4;
+        3
-    const ADAPTIVE: bool = true;
+    }
-    const DENSE: bool = true;
+
+    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
+            ode.y,  // y0 at start of step
-            next_y,          // y1 at end of step
+            next_y, // y1 at end of step
-            k[0],            // f(t0, y0) - derivative at start
+            k[0],   // f(t0, y0) - derivative at start
-            k[3],            // f(t1, y1) - derivative at end (FSAL)
+            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();
+        let bs3: BS3<1> = BS3::new();
-        assert_eq!(BS3::<1>::ORDER, 3);
+        assert_eq!(bs3.order(), 3);
-        assert_eq!(BS3::<1>::STAGES, 4);
+        assert_eq!(bs3.stages(), 4);
-        assert!(BS3::<1>::ADAPTIVE);
+        assert!(bs3.adaptive());
-        assert!(BS3::<1>::DENSE);
+        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.
+impl<const D: usize> Tolerances<D> for DormandPrince45<D>
-    pub fn a_tol(self, a_tol: f64) -> Self {
+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.
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
-    pub 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
+    fn r_tol(self, r_tol: f64) -> Self {
-    pub 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;
+    fn order(&self) -> usize {
-    const STAGES: usize = 7;
+        5
-    const ADAPTIVE: bool = true;
+    }
-    const DENSE: bool = true;
+
+    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,443 @@
+//! 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.
+    #[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
+    /// 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
+    /// `[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;
+    /// The order of the method's solution.
-    const STAGES: usize;
+    fn order(&self) -> usize;
-    const ADAPTIVE: bool;
-    const DENSE: bool;
 
-    /// Number of main stages stored in dense output (default: same as STAGES)
+    /// The number of stages the method evaluates per step.
-    const MAIN_STAGES: usize = Self::STAGES;
+    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)) / ε
+/// Uses forward differences: `J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε`
-/// where ε = √(machine_epsilon) * max(|y[j]|, 1.0)
+/// 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
+    /// Absolute tolerance for error estimation, per component
-    a_tol: f64,
+    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;
+    fn order(&self) -> usize {
-    const STAGES: usize = 2;
+        2
-    const ADAPTIVE: bool = true;
+    }
-    const DENSE: bool = true;
+
+    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 =
+            let scale = self.a_tol[i] + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs());
-                self.a_tol + 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 =
+            let rate = (errors[i].log10() - errors[i + 1].log10())
-                (errors[i].log10() - errors[i + 1].log10()) / (step_sizes[i].log10() - step_sizes[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);
+            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,
+        assert!(
-                "Convergence order {:.2} not within 0.2 of expected order 2", avg_order);
+            (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,
+        assert!(
-                "Adaptive solve should complete in less than 20 steps, got {}",
+            solution.states.len() < 20,
-                solution.states.len());
+            "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)
+impl<const D: usize> Tolerances<D> for Vern7<D>
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+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)
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
         self.a_tol = a_tol;
         self
     }
 
-    /// Set relative tolerance
+    fn r_tol(mut self, r_tol: f64) -> Self {
-    pub 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;
+    fn order(&self) -> usize {
-    const STAGES: usize = 10;
+        7
-    const ADAPTIVE: bool = true;
+    }
-    const DENSE: bool = true;
+
+    fn stages(&self) -> usize {
+        STAGES
+    }
+
+    fn adaptive(&self) -> bool {
+        true
+    }
+
+    fn dense(&self) -> bool {
+        true
+    }
 
     // Lazy dense output configuration
-    const MAIN_STAGES: usize = 10;
+    fn main_stages(&self) -> usize {
-    const EXTRA_STAGES: usize = 6;
+        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!(
+        assert!(ratio1 > 64.0 && ratio1 < 256.0, "First ratio: {}", ratio1);
-            ratio1 > 64.0 && ratio1 < 256.0,
+        assert!(ratio2 > 64.0 && ratio2 < 256.0, "Second ratio: {}", ratio2);
-            "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

@@ -1,8 +1,7 @@
 #![allow(dead_code)]
 // The readme is the crate's front page, and including it here means every
-// example in it is compiled and run by `cargo test`. Requires `autodiff`
+// example in it is compiled and run by `cargo test`.
-// because one of those examples differentiates.
+#![doc = include_str!("../readme.md")]
-#![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))]
 
 pub mod callback;
 pub mod controller;
@@ -16,10 +15,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 +51,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
@@ -159,7 +159,7 @@ mod tests {
 }
 
 /// Tests that the derivatives a dual instantiation produces are the right ones.
-#[cfg(all(test, feature = "autodiff"))]
+#[cfg(test)]
 mod autodiff_tests {
     use crate::ode::ProblemFunction;
     use crate::prelude::*;
@@ -321,6 +321,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 +375,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/scalar.rs

@@ -7,84 +7,33 @@
 //! differencing and no extra right-hand-side evaluations beyond the ones the
 //! dual arithmetic performs anyway.
 //!
-//! [`Scalar`] is implemented for `f64` always, and for the dual numbers of the
+//! [`Scalar`] extends the vocabulary [`num_dual`] already defines rather than
-//! [`num_dual`] crate under the `autodiff` feature (on by default).
+//! restating it: [`DualNum`] carries the arithmetic, the transcendentals and
+//! `From<f64>`, and [`DualStruct`] carries [`re`](DualStruct::re). Both are
+//! implemented for `f64` as well as for every dual type, so one blanket impl
+//! covers all of them, and a crate depending on this one alongside another
+//! built on `num-dual` can satisfy both with a single bound.
 
 use nalgebra::RealField;
+use num_dual::{DualNum, DualStruct};
 use num_traits::NumOps;
 
 /// A scalar that an ODE's state, time and step size may take.
-pub trait Scalar: RealField + Copy + NumOps<f64, Self> {
+///
-    /// This value with all derivative information discarded.
+/// `re`, which discards derivative information, is inherited from
-    fn re(&self) -> f64;
+/// [`DualStruct`]; this trait adds only [`lift`](Scalar::lift).
-
+pub trait Scalar:
+    RealField + Copy + NumOps<f64, Self> + DualNum<Primitive = f64> + DualStruct<Real = f64>
+{
     /// A real constant as a value of this type, with zero derivative.
-    fn lift(x: f64) -> Self;
-}
-
-impl Scalar for f64 {
-    #[inline(always)]
-    fn re(&self) -> f64 {
-        *self
-    }
-
-    #[inline(always)]
     fn lift(x: f64) -> Self {
-        x
+        Self::from(x)
     }
 }
 
-#[cfg(feature = "autodiff")]
+impl<T> Scalar for T where
-mod dual {
+    T: RealField + Copy + NumOps<f64, T> + DualNum<Primitive = f64> + DualStruct<Real = f64>
-    use super::Scalar;
+{
-    use nalgebra::{Const, RealField};
-    use num_dual::{DualNum, DualStruct};
-
-    /// Implements `Scalar` for one of `num-dual`'s dual number types.
-    macro_rules! impl_scalar {
-        ($ty:ident) => {
-            impl<T> Scalar for num_dual::$ty<T>
-            where
-                Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
-            {
-                #[inline(always)]
-                fn re(&self) -> f64 {
-                    DualStruct::re(self)
-                }
-
-                #[inline(always)]
-                fn lift(x: f64) -> Self {
-                    Self::from(x)
-                }
-            }
-        };
-        ($ty:ident, $($n:ident),+) => {
-            impl<T: nalgebra::Scalar, $(const $n: usize),+> Scalar
-                for num_dual::$ty<T, $(Const<$n>),+>
-            where
-                Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
-            {
-                #[inline(always)]
-                fn re(&self) -> f64 {
-                    DualStruct::re(self)
-                }
-
-                #[inline(always)]
-                fn lift(x: f64) -> Self {
-                    Self::from(x)
-                }
-            }
-        };
-    }
-
-    impl_scalar!(Dual);
-    impl_scalar!(Dual2);
-    impl_scalar!(Dual3);
-    impl_scalar!(HyperDual);
-    impl_scalar!(HyperHyperDual);
-    impl_scalar!(DualVec, N);
-    impl_scalar!(Dual2Vec, N);
-    impl_scalar!(HyperDualVec, M, N);
 }
 
 #[cfg(test)]
@@ -94,13 +43,12 @@ mod tests {
     #[test]
     fn test_f64_is_a_scalar() {
         assert_eq!(f64::lift(3.5), 3.5);
-        assert_eq!(Scalar::re(&3.5), 3.5);
+        assert_eq!(3.5f64.re(), 3.5);
     }
 
-    #[cfg(feature = "autodiff")]
     #[test]
     fn test_dual_projects_to_its_real_part() {
-        use num_dual::{Dual64, DualNum};
+        use num_dual::Dual64;
 
         // A lifted constant has no derivative; a seeded variable does, and both
         // project back to the same real value.
@@ -109,11 +57,10 @@ mod tests {
         assert_eq!(c.eps, 0.0);
 
         let x = Dual64::from_re(3.5).derivative();
-        assert_eq!(Scalar::re(&x), 3.5);
+        assert_eq!(x.re(), 3.5);
         assert_eq!(x.eps, 1.0);
 
         // d/dx x^3 = 3x^2 = 36.75 at x = 3.5
         assert_eq!(x.powi(3).eps, 3.0 * 3.5 * 3.5);
     }
-
 }