Merge branch 'release/0.8.0'

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().
2026-10-02 19:44:57 -04:00 · 2026-10-02 19:44:52 -04:00 · 2026-10-02 19:44:00 -04:00 · 2026-10-01 14:14:12 -04:00 · 2026-10-01 14:13:33 -04:00 · 2026-10-01 14:12:51 -04:00
29 changed files with 4085 additions and 821 deletions
@@ -1,25 +1,25 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.3.0"
+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://gitlab.rcjohnstone.com/connor/differential-equations"
+repository = "https://git.rcjohnstone.com/connor/ordinary-diffeq"
 license = "MIT"
-exclude = ["roadmap/"]
-
-# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+exclude = ["roadmap/", "hooks/"]

 [dependencies]
 serde = { version = "1.0", features = ["derive"] }
-nalgebra = { version = "0.34", features = ["serde-serialize"] }
+nalgebra = { version = "0.35", features = ["serde-serialize"] }
 num-traits = "0.2.19"
+num-dual = "0.15"
 roots = "0.0.8"

 [dev-dependencies]
 approx = "0.5"
 criterion = "0.7.0"
+serde_json = "1.0"

 [[bench]]
 name = "simple_1d"
@@ -1,17 +1,17 @@
 # Vern7 Performance Benchmark Report

-**Date**: 2025-10-24
-**Test System**: Linux 6.17.4-arch2-1
+**Date**: 2026-09-29
+**Test System**: Linux 7.2.6-arch2-1
 **Optimization Level**: Release build with full optimizations

 ## Executive Summary

 Vern7 demonstrates **substantial performance advantages** over lower-order methods (BS3 and DP5) at tight tolerances (1e-8 to 1e-12), achieving:
- **2.7x faster** than DP5 at 1e-10 tolerance (exponential problem)
- **3.8x faster** than DP5 in harmonic oscillator
- **8.8x faster** than DP5 for orbital mechanics
- **51x faster** than BS3 in harmonic oscillator
- **1.65x faster** than DP5 for interpolation workloads
+- **2.9x faster** than DP5 at 1e-10 tolerance (exponential problem)
+- **4.4x faster** than DP5 in harmonic oscillator
+- **1.1x faster** than DP5 for orbital mechanics
+- **83x faster** than BS3 in harmonic oscillator
+- **2.3x faster** than DP5 for interpolation workloads

 These results confirm Vern7's design goal: **maximum efficiency for high-accuracy requirements**.

@@ -23,14 +23,14 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac

 | Method | Time (μs) | Relative Speed | Speedup vs BS3 |
 |--------|-----------|----------------|----------------|
-| **Vern7** | **3.81** | **1.00x** (baseline) | **51.8x** |
-| DP5 | 10.43 | 2.74x slower | 18.9x |
-| BS3 | 197.37 | 51.8x slower | 1.0x |
+| **Vern7** | **2.96** | **1.00x** (baseline) | **53.2x** |
+| DP5 | 8.46 | 2.86x slower | 18.6x |
+| BS3 | 157.35 | 53.2x slower | 1.0x |

 **Analysis**:
- Vern7 is **2.7x faster** than DP5 and **51x faster** than BS3
+- Vern7 is **2.9x faster** than DP5 and **53x faster** than BS3
 - BS3's 3rd-order method requires many tiny steps to maintain 1e-10 accuracy
- DP5's 5th-order is better but still requires ~2.7x more work than Vern7
+- DP5's 5th-order is better but still requires ~2.9x more work than Vern7
 - Vern7's 7th-order allows much larger step sizes while maintaining accuracy

 ---
@@ -41,15 +41,15 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac

 | Method | Time (μs) | Relative Speed | Speedup vs BS3 |
 |--------|-----------|----------------|----------------|
-| **Vern7** | **26.89** | **1.00x** (baseline) | **55.1x** |
-| DP5 | 102.74 | 3.82x slower | 14.4x |
-| BS3 | 1,481.4 | 55.1x slower | 1.0x |
+| **Vern7** | **19.33** | **1.00x** (baseline) | **83.2x** |
+| DP5 | 84.19 | 4.36x slower | 19.1x |
+| BS3 | 1,608.7 | 83.2x slower | 1.0x |

 **Analysis**:
- Vern7 is **3.8x faster** than DP5 and **55x faster** than BS3
+- Vern7 is **4.4x faster** than DP5 and **83x faster** than BS3
 - Smooth periodic problems like harmonic oscillators are ideal for high-order methods
- BS3 requires ~1.5ms due to tiny steps needed for tight tolerance
- DP5 needs ~103μs, still significantly more than Vern7's 27μs
+- BS3 requires ~1.6ms due to tiny steps needed for tight tolerance
+- DP5 needs ~84μs, still significantly more than Vern7's 19μs
 - Higher dimensionality (2D vs 1D) amplifies the advantage of larger steps

 ---
@@ -60,15 +60,18 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac

 | Method | Time (μs) | Relative Speed | Speedup |
 |--------|-----------|----------------|---------|
-| **Vern7** | **98.75** | **1.00x** (baseline) | **8.77x** |
-| DP5 | 865.79 | 8.77x slower | 1.0x |
+| **Vern7** | **83.88** | **1.00x** (baseline) | **1.13x** |
+| DP5 | 94.68 | 1.13x slower | 1.0x |

 **Analysis**:
- Vern7 is **8.8x faster** than DP5 for this challenging 6D problem
+- Vern7 is **1.1x faster** than DP5 on this problem
+- Step counts: Vern7 184, DP5 440. Vern7 takes 2.4x fewer steps, but each step
+  costs 16 stages against DP5's 7, so the advantage very nearly cancels
 - Orbital mechanics requires tight tolerances to maintain energy conservation
 - BS3 was too slow to include in the benchmark at this tolerance
- 6D problem with long integration time shows Vern7's scalability
- This represents realistic astrodynamics/orbital mechanics workloads
+- Vern7's advantage here is in step count, which is what matters once the
+  right-hand side is expensive - a real gravity model rather than a two-body
+  point mass

 ---

@@ -78,17 +81,20 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac

 | Method | Time (μs) | Relative Speed | Notes |
 |--------|-----------|----------------|-------|
-| **Vern7** | **11.05** | **1.00x** (baseline) | Lazy extra stages |
-| DP5 | 18.27 | 1.65x slower | Standard dense output |
+| **Vern7** | **7.28** | **1.00x** (baseline) | Lazy extra stages |
+| DP5 | 16.73 | 2.30x slower | Standard dense output |

 **Analysis**:
- Vern7 with lazy computation is **1.65x faster** than DP5
+- Vern7 with lazy computation is **2.3x faster** than DP5
 - First interpolation triggers lazy computation of 6 extra stages (k11-k16)
 - Subsequent interpolations reuse cached extra stages (~10ns RefCell overhead)
 - Despite computing extra stages, Vern7 is still faster overall due to:
  1. Fewer total integration steps (larger step sizes)
  2. Higher accuracy interpolation (7th order vs 5th order)
- Lazy computation adds minimal overhead (~6μs for 6 stages, amortized over 100 interpolations)
+- Lazy computation adds minimal overhead, amortized over 100 interpolations
+- Each benchmark iteration builds a fresh `Solution`, so the lazy stages are
+  computed once per iteration rather than served from the `RefCell` cache - the
+  cost of interpolating into a solution you have just produced

 ---

@@ -100,11 +106,11 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac

 | Tolerance | DP5 (μs) | Vern7 (μs) | Speedup | Winner |
 |-----------|----------|------------|---------|--------|
-| 1e-6 | 2.63 | 2.05 | 1.28x | Vern7 |
-| 1e-7 | 3.71 | 2.74 | 1.35x | Vern7 |
-| 1e-8 | 5.43 | 3.12 | 1.74x | Vern7 |
-| 1e-9 | 7.97 | 3.86 | 2.06x | **Vern7** |
-| 1e-10 | 11.33 | 5.33 | 2.13x | **Vern7** |
+| 1e-6 | 1.94 | 1.51 | 1.29x | Vern7 |
+| 1e-7 | 2.82 | 1.93 | 1.46x | Vern7 |
+| 1e-8 | 4.20 | 2.29 | 1.83x | Vern7 |
+| 1e-9 | 6.25 | 2.97 | 2.11x | **Vern7** |
+| 1e-10 | 9.52 | 3.89 | 2.45x | **Vern7** |

 ### Performance Scaling Chart (Conceptual)

@@ -128,7 +134,7 @@ Time (μs)

 **Analysis**:
 - At **moderate tolerances (1e-6)**: Vern7 is 1.3x faster
- At **tight tolerances (1e-10)**: Vern7 is 2.1x faster
+- At **tight tolerances (1e-10)**: Vern7 is 2.5x faster
 - **Crossover point**: Vern7 becomes increasingly advantageous as tolerance tightens
 - DP5's time scales roughly quadratically with tolerance
 - Vern7's time scales more slowly (higher order = larger steps)
@@ -157,22 +163,23 @@ Time (μs)
 ### Lazy Computation Impact

 The lazy computation of extra stages (k11-k16) provides:
- **Minimal overhead**: ~6μs to compute 6 extra stages
+- **Minimal overhead**: the 6 extra stages are computed once per step that is interpolated into
 - **Cache efficiency**: Extra stages computed once per interval, reused for multiple interpolations
 - **Memory efficiency**: Only computed when interpolation is requested
- **Performance**: Despite extra computation, still 1.65x faster than DP5 for interpolation workloads
+- **Performance**: Despite extra computation, still 2.3x faster than DP5 for interpolation workloads

 ### Step Size Comparison

-Estimated step sizes at 1e-10 tolerance for exponential problem:
+Measured at 1e-10 tolerance on the exponential problem (`y' = y`, t = 0 to 4):

 | Method | Avg Step Size | Steps Required | Function Evals |
 |--------|---------------|----------------|----------------|
-| BS3 | ~0.002 | ~2000 | ~8000 |
-| DP5 | ~0.01 | ~400 | ~2400 |
-| **Vern7** | ~0.05 | **~80** | **~800** |
+| BS3 | 0.0014 | 2,916 | 11,664 |
+| DP5 | 0.0354 | 113 | 791 |
+| **Vern7** | 0.1379 | **29** | **290** |

-**Vern7 requires ~3x fewer function evaluations than DP5.**
+**Vern7 requires 2.7x fewer function evaluations than DP5**, taking 3.9x fewer
+steps at 10 stages each against DP5's 7.

 ---

@@ -185,7 +192,7 @@ Our Rust implementation achieves performance comparable to Julia's highly-optimi
 | Step computation | Highly optimized, FSAL | Optimized, no FSAL |
 | Lazy interpolation | ✓ | ✓ |
 | Stage caching | RefCell-based | RefCell-based (~10ns) |
-| Memory allocation | Minimal | Minimal |
+| Memory allocation | Cached in the integrator | Two heap allocations per step |
 | Relative speed | Baseline | ~Comparable |

 **Note**: Direct comparison difficult due to different hardware and problems, but algorithmic approach is identical.
@@ -213,8 +220,8 @@ Our Rust implementation achieves performance comparable to Julia's highly-optimi

 Vern7 successfully achieves its design goal of being the **most efficient method for high-accuracy non-stiff problems**. The implementation with lazy computation of extra stages provides:

- ✅ **2-9x speedup** over DP5 at tight tolerances
- ✅ **50x+ speedup** over BS3 at tight tolerances
+- ✅ **1.1-4.4x speedup** over DP5 at tight tolerances
+- ✅ **53-83x speedup** over BS3 at tight tolerances
 - ✅ **Efficient lazy interpolation** with minimal overhead
 - ✅ **Full 7th-order accuracy** for both steps and interpolation
 - ✅ **Memory-efficient caching** with RefCell
@@ -225,7 +232,7 @@ The results validate the effort invested in implementing the complex 16-stage in

 ## Appendix: Benchmark Configuration

-**Hardware**: Not specified (Linux system)
+**Hardware**: Linux 7.2.6-arch2-1
 **Compiler**: rustc (release mode, full optimizations)
 **Measurement Tool**: Criterion.rs v0.7.0
 **Sample Size**: 100 samples per benchmark
@@ -1,6 +1,6 @@
 # BS3 vs DP5 Benchmark Results

-Generated: 2025-10-23
+Generated: 2026-09-29

 ## Summary

@@ -10,7 +10,7 @@ Comprehensive performance comparison between **BS3** (Bogacki-Shampine 3rd order

 ### Overall Performance Comparison

-**DP5 is consistently faster than BS3 across all tested scenarios**, typically by a factor of **1.5x to 4.3x**.
+**DP5 is consistently faster than BS3 across all tested scenarios**, typically by a factor of **1.2x to 5.1x**.

 This might seem counterintuitive since BS3 uses fewer stages (4 vs 7), but several factors explain DP5's superior performance:

@@ -32,8 +32,8 @@ Despite being slower in these benchmarks, BS3 still has value:

 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 3.28 µs | 1.92x slower |
-| **DP5** | 1.70 µs | baseline |
+| **BS3** | 2.99 µs | 2.21x slower |
+| **DP5** | 1.35 µs | baseline |

 Simple 1D problem with smooth exponential solution.

@@ -41,8 +41,8 @@ Simple 1D problem with smooth exponential solution.

 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 30.70 µs | 2.25x slower |
-| **DP5** | 13.67 µs | baseline |
+| **BS3** | 27.37 µs | 2.76x slower |
+| **DP5** | 9.92 µs | baseline |

 2D conservative system with periodic solution.

@@ -50,8 +50,8 @@ Simple 1D problem with smooth exponential solution.

 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 132.35 µs | 3.57x slower |
-| **DP5** | 37.11 µs | baseline |
+| **BS3** | 113.44 µs | 3.81x slower |
+| **DP5** | 29.74 µs | baseline |

 Nonlinear 2D system with trigonometric terms.

@@ -59,8 +59,8 @@ Nonlinear 2D system with trigonometric terms.

 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 124.72 µs | 1.45x slower |
-| **DP5** | 86.10 µs | baseline |
+| **BS3** | 98.02 µs | 5.07x slower |
+| **DP5** | 19.32 µs | baseline |

 Higher-dimensional problem with gravitational dynamics.

@@ -68,8 +68,8 @@ Higher-dimensional problem with gravitational dynamics.

 | Method | Time (solve + 100 interpolations) | Ratio |
 |--------|-----------------------------------|-------|
-| **BS3** | 19.68 µs | 4.81x slower |
-| **DP5** | 4.09 µs | baseline |
+| **BS3** | 16.34 µs | 4.67x slower |
+| **DP5** | 3.50 µs | baseline |

 BS3 uses cubic Hermite interpolation, DP5 uses optimized 5th order interpolation.

@@ -79,14 +79,11 @@ Performance across different tolerance levels (`y' = -y` problem):

 | Tolerance | BS3 Time | DP5 Time | Ratio (BS3/DP5) |
 |-----------|----------|----------|-----------------|
-| 1e-3 | 1.63 µs | 1.26 µs | 1.30x |
-| 1e-4 | 2.61 µs | 1.54 µs | 1.70x |
-| 1e-5 | 4.64 µs | 2.03 µs | 2.28x |
-| 1e-6 | 8.76 µs | ~2.6 µs* | ~3.4x* |
-| 1e-7 | -** | -** | - |
-
-\* Estimated from trend (benchmark timed out)
-\** Not completed
+| 1e-3 | 1.10 µs | 0.91 µs | 1.21x |
+| 1e-4 | 1.87 µs | 1.19 µs | 1.57x |
+| 1e-5 | 3.46 µs | 1.44 µs | 2.40x |
+| 1e-6 | 6.77 µs | 1.95 µs | 3.46x |
+| 1e-7 | 13.73 µs | 2.79 µs | 4.91x |

 **Observation**: The performance gap widens as tolerance tightens, because DP5's higher order allows it to take larger steps while maintaining accuracy.

@@ -96,15 +93,18 @@ Performance across different tolerance levels (`y' = -y` problem):

 1. **DP5 is the better default choice** for most problems requiring moderate to high accuracy
 2. **Performance gap increases** with tighter tolerances (favoring DP5)
-3. **Higher dimensions** slightly favor BS3 relative to DP5 (1.45x vs 3.57x slowdown)
-4. **Interpolation** strongly favors DP5 (4.8x faster)
+3. **Higher dimensions favor DP5 most strongly** (5.07x on the 6D orbit against
+   2.21x on 1D exponential decay)
+4. **Interpolation** strongly favors DP5 (4.7x faster)

 ### Implementation Quality

 Both integrators pass all accuracy and convergence tests:
 - ✅ BS3: 3rd order convergence rate verified
 - ✅ DP5: 5th order convergence rate verified (existing implementation)
- ✅ Both: FSAL property correctly implemented
+- ⚠️ Both: FSAL is a property of the *methods*, but neither implementation
+  exploits it - each step recomputes `k[0] = f(t, y)` rather than reusing the
+  previous step's last stage. See roadmap feature 36.
 - ✅ Both: Dense output accurate to specified order

 ### Future Optimizations
@@ -27,9 +27,8 @@ fn bench_exponential_decay(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10.0, y0, params);
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, bs3, controller).solve();
-            });
+            let mut problem = Problem::new(ode, bs3, controller);
+            black_box(problem.solve());
        });
    });

@@ -37,9 +36,8 @@ fn bench_exponential_decay(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10.0, y0, params);
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -66,9 +64,8 @@ fn bench_harmonic_oscillator(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 20.0, y0, ());
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, bs3, controller).solve();
-            });
+            let mut problem = Problem::new(ode, bs3, controller);
+            black_box(problem.solve());
        });
    });

@@ -76,9 +73,8 @@ fn bench_harmonic_oscillator(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 20.0, y0, ());
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -109,9 +105,8 @@ fn bench_pendulum(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10.0, y0, params);
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, bs3, controller).solve();
-            });
+            let mut problem = Problem::new(ode, bs3, controller);
+            black_box(problem.solve());
        });
    });

@@ -119,9 +114,8 @@ fn bench_pendulum(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10.0, y0, params);
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -160,9 +154,8 @@ fn bench_orbit_6d(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params);
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, bs3, controller).solve();
-            });
+            let mut problem = Problem::new(ode, bs3, controller);
+            black_box(problem.solve());
        });
    });

@@ -170,9 +163,8 @@ fn bench_orbit_6d(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params);
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -199,11 +191,14 @@ fn bench_interpolation(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                let solution = Problem::new(ode, bs3, controller).solve();
-                // Interpolate at 100 points
-                let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect();
-            });
+            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<_>>(),
+            );
        });
    });

@@ -212,11 +207,14 @@ fn bench_interpolation(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                let solution = Problem::new(ode, dp45, controller).solve();
-                // Interpolate at 100 points
-                let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect();
-            });
+            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<_>>(),
+            );
        });
    });

@@ -243,9 +241,8 @@ fn bench_tolerance_scaling(c: &mut Criterion) {
            let ode = ODE::new(&derivative, 0.0, 10.0, y0, ());
            let bs3 = BS3::new().a_tol(tol).r_tol(tol);
            b.iter(|| {
-                black_box({
-                    Problem::new(ode, bs3, controller).solve();
-                });
+                let mut problem = Problem::new(ode, bs3, controller);
+                black_box(problem.solve());
            });
        });

@@ -253,9 +250,8 @@ fn bench_tolerance_scaling(c: &mut Criterion) {
            let ode = ODE::new(&derivative, 0.0, 10.0, y0, ());
            let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
            b.iter(|| {
-                black_box({
-                    Problem::new(ode, dp45, controller).solve();
-                });
+                let mut problem = Problem::new(ode, dp45, controller);
+                black_box(problem.solve());
            });
        });
    }
@@ -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;
@@ -29,9 +29,8 @@ fn bench_orbit(c: &mut Criterion) {

    c.bench_function("bench_orbit", |b| {
        b.iter(|| {
-            std::hint::black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            std::hint::black_box(problem.solve());
        });
    });
 }
@@ -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,);
@@ -20,9 +20,8 @@ fn bench_simple_1d(c: &mut Criterion) {

    c.bench_function("bench_simple_1d", |b| {
        b.iter(|| {
-            std::hint::black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            std::hint::black_box(problem.solve());
        });
    });
 }
@@ -44,10 +43,13 @@ fn bench_interpolation_1d(c: &mut Criterion) {

    c.bench_function("bench_interpolation_1d", |b| {
        b.iter(|| {
-            std::hint::black_box({
-                let solution = Problem::new(ode, dp45, controller).solve();
-                let _ = (0..100).map(|t| solution.interpolate(t as f64 * 0.1)[0]);
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            let solution = problem.solve();
+            std::hint::black_box(
+                (0..100)
+                    .map(|t| solution.interpolate(t as f64 * 0.1)[0])
+                    .sum::<f64>(),
+            );
        });
    });
 }
@@ -28,9 +28,8 @@ fn bench_exponential_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 4.0, y0, ());
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, bs3, controller).solve();
-            });
+            let mut problem = Problem::new(ode, bs3, controller);
+            black_box(problem.solve());
        });
    });

@@ -38,9 +37,8 @@ fn bench_exponential_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 4.0, y0, ());
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -48,9 +46,8 @@ fn bench_exponential_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 4.0, y0, ());
        let vern7 = Vern7::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, vern7, controller).solve();
-            });
+            let mut problem = Problem::new(ode, vern7, controller);
+            black_box(problem.solve());
        });
    });

@@ -77,9 +74,8 @@ fn bench_harmonic_oscillator_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 20.0, y0, ());
        let bs3 = BS3::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, bs3, controller).solve();
-            });
+            let mut problem = Problem::new(ode, bs3, controller);
+            black_box(problem.solve());
        });
    });

@@ -87,9 +83,8 @@ fn bench_harmonic_oscillator_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 20.0, y0, ());
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -97,9 +92,8 @@ fn bench_harmonic_oscillator_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 20.0, y0, ());
        let vern7 = Vern7::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, vern7, controller).solve();
-            });
+            let mut problem = Problem::new(ode, vern7, controller);
+            black_box(problem.solve());
        });
    });

@@ -138,9 +132,8 @@ fn bench_orbit_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params);
        let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, dp45, controller).solve();
-            });
+            let mut problem = Problem::new(ode, dp45, controller);
+            black_box(problem.solve());
        });
    });

@@ -148,9 +141,8 @@ fn bench_orbit_tight_tol(c: &mut Criterion) {
        let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params);
        let vern7 = Vern7::new().a_tol(tol).r_tol(tol);
        b.iter(|| {
-            black_box({
-                Problem::new(ode, vern7, controller).solve();
-            });
+            let mut problem = Problem::new(ode, vern7, controller);
+            black_box(problem.solve());
        });
    });

@@ -175,27 +167,31 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
    // Vern7 with interpolation (should compute extra stages lazily)
    group.bench_function("vern7_with_interpolation", |b| {
        b.iter(|| {
-            black_box({
-                let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
-                let vern7 = Vern7::new().a_tol(tol).r_tol(tol);
-                let mut problem = Problem::new(ode, vern7, controller);
-                let solution = problem.solve();
-                // Interpolate at 100 points - first one computes extra stages
-                let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect();
-            });
+            let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
+            let vern7 = Vern7::new().a_tol(tol).r_tol(tol);
+            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<_>>(),
+            );
        });
    });

    // DP5 with interpolation for comparison
    group.bench_function("dp5_with_interpolation", |b| {
        b.iter(|| {
-            black_box({
-                let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
-                let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
-                let mut problem = Problem::new(ode, dp45, controller);
-                let solution = problem.solve();
-                let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect();
-            });
+            let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
+            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<_>>(),
+            );
        });
    });

@@ -223,9 +219,8 @@ fn bench_tolerance_scaling_vern7(c: &mut Criterion) {
            let ode = ODE::new(&derivative, 0.0, 10.0, y0, ());
            let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
            b.iter(|| {
-                black_box({
-                    Problem::new(ode, dp45, controller).solve();
-                });
+                let mut problem = Problem::new(ode, dp45, controller);
+                black_box(problem.solve());
            });
        });

@@ -233,9 +228,8 @@ fn bench_tolerance_scaling_vern7(c: &mut Criterion) {
            let ode = ODE::new(&derivative, 0.0, 10.0, y0, ());
            let vern7 = Vern7::new().a_tol(tol).r_tol(tol);
            b.iter(|| {
-                black_box({
-                    Problem::new(ode, vern7, controller).solve();
-                });
+                let mut problem = Problem::new(ode, vern7, controller);
+                black_box(problem.solve());
            });
        });
    }
@@ -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
@@ -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
@@ -1,4 +1,4 @@
-# DifferentialEquations
+# ordinary-diffeq

 A library, written in Rust, for integrating ordinary differential equations. For now, this is
 relatively simple, but it does have key features that are needed for orbit propagation, ray tracing,
@@ -15,18 +15,22 @@ and field line tracing:
 | **Vern7** (Verner) | 7(6) | 10+6 | 7th order | High accuracy (rtol ~ 1e-8 to 1e-12) |

 **Performance at 1e-10 tolerance:**
- Vern7: **2.7-8.8x faster** than DP5
- Vern7: **50x+ faster** than BS3
+- Vern7: **1.1-4.4x faster** than DP5
+- Vern7: **53-83x faster** than BS3

 See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analysis.

 ### 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

@@ -47,7 +51,7 @@ pendulum problem):

 ```rust
 use nalgebra::Vector2;
-use differential_equations::prelude::*;
+use ordinary_diffeq::prelude::*;
 use std::f64::consts::PI;

 // Define the system (parameters, derivative, and initial state)
@@ -68,15 +72,211 @@ 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 = Callback {
-    event: &|t: f64, _y: Vector2<f64>, _p: &Params| { 5.0 - t },
+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, which is a blanket extension of
+`DualNum` and `DualStruct` from [`num-dual`](https://crates.io/crates/num-dual): those two supply
+the arithmetic, the transcendentals and `re`, which discards derivative information, and `Scalar`
+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:
+
+```rust,ignore
+fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
+    let &(g, l) = p;
+    Vector2::new( y[1], -(g/l) * y[0].sin() )
+}
+
+fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
+    let &(g, l) = p;
+    Vector2::new( y[1], -y[0].sin() * (g/l) )
+}
+```
+
+`Params` stays `(f64, f64)` — parameters only need to become `T` themselves if you want to
+differentiate *with respect to* them. Real constants mix into `T` arithmetic directly; `Scalar::lift`
+is only needed to scale a vector by a real, or to build a `T` out of nothing.
+
+With that one change, the same function can be solved or differentiated:
+
+```rust
+use nalgebra::{Matrix2, Vector2};
+use num_dual::DualSVec64;
+use ordinary_diffeq::prelude::*;
+use std::f64::consts::PI;
+
+type Params = (f64, f64);
+
+fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
+    let &(g, l) = p;
+    Vector2::new(y[1], -y[0].sin() * (g / l))
+}
+
+// One propagation in `DualSVec64<2>` yields the whole 2x2 Jacobian of the
+// final state with respect to the initial one.
+let (y_end, jacobian): (Vector2<f64>, Matrix2<f64>) = num_dual::jacobian(
+    |y0: Vector2<DualSVec64<2>>| {
+        type T = DualSVec64<2>;
+        let ode = ODE::new(&derivative, T::lift(0.0), T::lift(1.0), y0, (9.81, 1.0));
+        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();
+        *solution.states.last().unwrap()
+    },
+    &Vector2::new(0.0, PI / 2.0),
+);
+
+// The pendulum is Hamiltonian, so its flow is area-preserving
+assert!((jacobian.determinant() - 1.0).abs() < 1e-9);
+```
+
+These are exact derivatives of the scheme as it ran, at machine precision, with no finite
+differencing and no step size to tune.
@@ -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:
@@ -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
@@ -298,6 +332,12 @@ Each feature below links to a detailed implementation plan in the `features/` di
  - **Dependencies**: None
  - **Effort**: Small

+- [ ] **[Reduce Per-Step Allocation](features/39-per-step-allocation.md)**
+  - Two heap allocations per accepted step (stage array, dense coefficients)
+  - Against ~150 ns of arithmetic; matters only for cheap right-hand sides
+  - **Dependencies**: None
+  - **Effort**: Medium
+
 ---

 ## Implementation Notes
@@ -330,17 +370,20 @@ Each algorithm implementation should include:

 ## Progress Tracking

-Total Features: 38
- Tier 1: 8 features (4/8 complete) ✅
- Tier 2: 12 features (0/12 complete)
- Tier 3: 18 features (0/18 complete)
+Total Features: 42
+- Tier 1: 11 features (7/11 complete) ✅
+- Tier 2: 11 features (0/11 complete)
+- Tier 3: 20 features (0/20 complete)

-**Overall Progress: 10.5% (4/38 features complete)**
+**Overall Progress: 16.7% (7/42 features complete)**

 ### Completed Features
 1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23)
 2. ✅ Vern7 (Verner 7th order) - Tier 1 (2025-10-24)
 3. ✅ PID Controller - Tier 1 (2025-10-24)
 4. ✅ Rosenbrock23 - Tier 1 (2025-10-24)
+5. ✅ Tolerances Trait - Tier 1 (2026-09-30)
+6. ✅ Runtime Integrator Selection - Tier 1 (2026-09-30)
+7. ✅ Callback Effects That Can Change the State - Tier 1 (2026-10-01)

-Last updated: 2025-10-24
+Last updated: 2026-10-01
@@ -4,6 +4,31 @@

 Finite difference and auto-diff Jacobians

+## Note - 2026-09-29
+
+Forward-mode AD arrived by a different route than this entry assumes. On the
+`feature/generic-scalar` branch, state, time and step size became generic over a
+`Scalar` trait, implemented by `f64` and by `num-dual`'s dual numbers. That gives
+differentiation *through the whole solve* - d y(t_f)/d y0, d/d t0, d/d t_end and
+d/d params, so a full 6x6 state transition matrix comes out of one propagation.
+This entry is framed around the narrower case, df/dy at a point for implicit
+methods, which now falls out of the same machinery.
+
+What that leaves:
+
+- `rosenbrock::compute_jacobian` is still forward finite differences
+  (eps = sqrt(machine_eps) * max(|y[j]|, 1), D+1 evaluations, ~8 digits).
+  Pointing it at a dual evaluation is now a small change, not a large one.
+- Sparse Jacobian support: untouched.
+- The **Effort: Large** estimate below predates the genericisation and wants
+  revisiting.
+
+Deferred deliberately, as scope creep against the genericisation: how a Jacobian
+should be obtained for whichever algorithms need one. Whether the integrator
+requests it through a trait, whether AD or finite differences is a per-problem
+choice, and how either interacts with the `JacobianUpdateStrategy` reuse policy
+already in `rosenbrock.rs`. To be thought through properly on its own.
+
 ## Why This Feature Matters

 (To be detailed)
@@ -0,0 +1,81 @@
+# Feature: Reduce Per-Step Allocation
+
+## Overview
+
+Remove the heap allocations every accepted step performs.
+
+## Why This Feature Matters
+
+Each call to `step()` allocates a `Vec` for the stage values
+(`vec![SVector::zeros(); STAGES]`) and a second `Vec` for the dense output
+coefficients. `Problem::solve` then wraps the latter in a `RefCell` and pushes
+it onto a growing vector. That is two allocations per step against roughly
+150 ns of arithmetic.
+
+Measured 2026-09-29 on the harmonic oscillator at 1e-10, where the right-hand
+side is two flops so the timing is almost entirely integrator overhead:
+
+| Method | Steps | Total | Per step |
+|--------|-------|-------|----------|
+| BS3 | 13,927 | 1.58 ms | 113 ns |
+| DP45 | 722 | 112 us | 155 ns |
+| Vern7 | 124 | 26.6 us | 215 ns |
+
+`BS3::step` in isolation is 44.6 ns for four stages.
+
+The gain is largest exactly where it is least important: problems with cheap
+right-hand sides. For the orbit propagation this crate is aimed at, a
+spherical-harmonic gravity evaluation dominates and per-step overhead is
+rounding error. Worth doing, not worth doing first.
+
+## Dependencies
+
+None. `STAGES` is already an associated const, so the stage array can be a
+fixed-size array rather than a `Vec`.
+
+## Implementation Approach
+
+(To be detailed.) Sketch: stage storage becomes `[SVector<T, D>; STAGES]` if
+const generics allow it at the trait level, or a caller-supplied scratch buffer
+threaded through `step()`. Dense coefficients want a pre-reserved buffer sized
+from an estimate of the step count, or a flat `Vec<SVector<T, D>>` indexed by
+stride rather than a `Vec` of `Vec`s.
+
+Note the interaction with the `RefCell` in `Solution::dense`: it exists so
+Vern7 can append its lazily computed extra stages. A flat buffer has to keep
+that possible.
+
+## Implementation Tasks
+
+- [ ] Core implementation
+- [ ] Integration with existing code
+- [ ] Testing
+- [ ] Benchmarking
+
+## Testing Requirements
+
+Numerics must not change. The existing benches plus a step-count and
+final-state comparison against the previous release are sufficient: results
+should be bit-identical.
+
+## References
+
+1. Julia implementation: OrdinaryDiffEq.jl, which caches stage storage in the
+   integrator rather than allocating per step
+
+## Complexity Estimate
+
+**Effort**: Medium
+
+**Risk**: Low - no change to the numerics, so the test suite is a strong guard
+
+## Success Criteria
+
+- [ ] No allocation in the steady-state step loop
+- [ ] Bit-identical results against the previous release
+- [ ] Measurable improvement on the cheap-RHS benchmarks
+
+## Future Enhancements
+
+Supersedes the "Optimize memory allocation in step method" bullet under Future
+Enhancements in `01-bs3-method.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
@@ -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
@@ -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
@@ -1,22 +1,120 @@
+//! Conditions that interrupt an integration, and what they do when they fire.
+
 use nalgebra::SVector;

-use super::ode::ODE;
+use crate::ode::ODE;
+use crate::scalar::Scalar;

-/// A function that takes in a time and a state and outputs a single float value
+/// A condition the solver watches for, and the change it makes when it is met.
 ///
-/// The integration solver will check this function for zero crossings
-#[derive(Clone, Copy)]
-pub struct Callback<'a, const D: usize, P> {
-    /// The function to check for zero crossings
-    pub event: &'a dyn Fn(f64, SVector<f64, D>, &P) -> f64,
+/// The solver evaluates [`event`](Callback::event) at each end of every accepted
+/// step and looks for a sign change. When it finds one it shortens the step to
+/// land on the crossing, records the state there, and then calls
+/// [`effect`](Callback::effect) on the ODE. Whatever the effect leaves in
+/// `ode.y` is the state the integration carries on from.
+///
+/// An effect may:
+///
+/// - change the state, which records a jump: the solution holds both the value
+///   at the crossing and the value the effect left, at the same time (see
+///   [`Solution`](crate::problem::Solution));
+/// - change the parameters, which changes the right-hand side from the crossing
+///   onward;
+/// - end the solution at the crossing, by calling [`stop`] or
+///   [`ODE::terminate`];
+/// - record something about itself, since it takes `&mut self`.
+///
+/// # Example
+///
+/// A fire-once callback, which needs the `&mut self` on the effect to remember
+/// that it has run:
+///
+/// ```
+/// use nalgebra::Vector1;
+/// use ordinary_diffeq::prelude::*;
+///
+/// struct Kick {
+///     fired: usize,
+/// }
+///
+/// impl Callback<f64, 1, ()> for Kick {
+///     fn event(&self, _t: f64, y: &Vector1<f64>, _p: &()) -> f64 {
+///         5.0 - y[0]
+///     }
+///
+///     fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) {
+///         self.fired += 1;
+///         ode.y = Vector1::new(100.0);
+///     }
+/// }
+///
+/// fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+///     Vector1::new(1.0)
+/// }
+///
+/// let ode = ODE::new(&derivative, 0.0, 10.0, Vector1::new(0.0), ());
+/// let mut kick = Kick { fired: 0 };
+/// let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+///     .with_callback(&mut kick);
+/// let solution = problem.solve();
+///
+/// // The state jumped to 100 at t = 5 and kept integrating from there.
+/// assert!((solution.states.last().unwrap()[0] - 105.0).abs() < 1e-9);
+/// ```
+pub trait Callback<T: Scalar, const D: usize, P> {
+    /// The function whose zero crossings mark the event.
+    ///
+    /// Called repeatedly while the crossing is being located, so it must not
+    /// depend on how many times it has been called.
+    fn event(&self, t: T, y: &SVector<T, D>, params: &P) -> T;

-    /// The function to change the ODE
-    pub effect: &'a dyn Fn(&mut ODE<D, P>),
+    /// 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 convenience function for stopping the integration
-pub fn stop<const D: usize, P>(ode: &mut ODE<D, P>) {
-    ode.t_end = ode.t;
+/// 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 mut too_high = FnCallback {
+///     event: &|_t: f64, y: &Vector2<f64>, _p: &Params| 10.0 - y[0],
+///     effect: &stop,
+/// };
+/// ```
+#[derive(Clone, Copy)]
+pub struct FnCallback<'a, T: Scalar, const D: usize, P> {
+    /// The function to check for zero crossings
+    pub event: &'a dyn Fn(T, &SVector<T, D>, &P) -> T,
+
+    /// The function to change the ODE
+    pub effect: &'a dyn Fn(&mut ODE<'_, T, D, P>),
+}
+
+impl<T: Scalar, const D: usize, P> Callback<T, D, P> for FnCallback<'_, T, D, P> {
+    #[inline]
+    fn event(&self, t: T, y: &SVector<T, D>, params: &P) -> T {
+        (self.event)(t, y, params)
+    }
+
+    #[inline]
+    fn effect(&mut self, ode: &mut ODE<'_, T, D, P>) {
+        (self.effect)(ode)
+    }
+}
+
+/// 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.terminate();
 }

 #[cfg(test)]
@@ -26,8 +124,8 @@ mod tests {
    #[test]
    fn test_basic_callbacks() {
        type Params = ();
-        let _value_too_high = Callback {
-            event: &|_: f64, y: SVector<f64, 3>, _p: &Params| 10.0 - y[0],
+        let _value_too_high = FnCallback {
+            event: &|_: f64, y: &SVector<f64, 3>, _p: &Params| 10.0 - y[0],
            effect: &stop,
        };
    }
@@ -1,11 +1,13 @@
+use crate::scalar::Scalar;
+
 #[derive(Debug, Clone, Copy, PartialEq)]
-pub enum TryStep {
-    Accepted(f64, f64),
-    NotYetAccepted(f64),
+pub enum TryStep<T> {
+    Accepted(T, T),
+    NotYetAccepted(T),
 }

-impl TryStep {
-    pub fn extract(&self) -> f64 {
+impl<T: Copy> TryStep<T> {
+    pub fn extract(&self) -> T {
        match self {
            TryStep::Accepted(h, _) => *h,
            TryStep::NotYetAccepted(h) => *h,
@@ -16,7 +18,7 @@ impl TryStep {
        matches!(self, TryStep::Accepted(_, _))
    }

-    pub fn reset(&mut self) -> Result<TryStep, &str> {
+    pub fn reset(&mut self) -> Result<TryStep<T>, &str> {
        match self {
            TryStep::Accepted(_, h) => Ok(TryStep::NotYetAccepted(*h)),
            TryStep::NotYetAccepted(_) => Err("Cannot reset a NotYetAccepted TryStep"),
@@ -24,12 +26,26 @@ impl TryStep {
    }
 }

-pub trait Controller<const D: usize> {
-    fn determine_step(&mut self, h: f64, err: f64) -> TryStep;
+/// Chooses the next step size from the last one and its normalized error.
+pub trait Controller<T: Scalar> {
+    fn determine_step(&mut self, h: T, err: f64) -> TryStep<T>;
+
+    /// The step size to attempt next
+    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)]
-pub struct PIController {
+pub struct PIController<T: Scalar = f64> {
    pub alpha: f64,
    pub beta: f64,
    pub factor_c1: f64,
@@ -37,13 +53,13 @@ pub struct PIController {
    pub factor_old: f64,
    pub h_max: f64,
    pub safety_factor: f64,
-    pub next_step_guess: TryStep,
+    pub next_step_guess: TryStep<T>,
 }

-impl<const D: usize> Controller<D> for PIController {
+impl<T: Scalar> Controller<T> for PIController<T> {
    /// Determines if the previously run step size and error were valid or not. Either way, it also
    /// returns what the next step size should be
-    fn determine_step(&mut self, prev_step: f64, err: f64) -> TryStep {
+    fn determine_step(&mut self, prev_step: T, err: f64) -> TryStep<T> {
        let factor_11 = err.powf(self.alpha);
        let factor = self.factor_c2.max(
            self.factor_c1
@@ -53,19 +69,31 @@ impl<const D: usize> Controller<D> for PIController {
            let mut h = prev_step / factor;
            // Accept the stepsize and provide what the next step size should be
            self.factor_old = err.max(1.0e-4);
-            if h.abs() > self.h_max {
+            if h.re().abs() > self.h_max {
                // If the step goes past the maximum allowed, though, we shrink it
-                h = self.h_max.copysign(h);
+                h = T::lift(self.h_max.copysign(h.re()));
            }
            TryStep::Accepted(prev_step, h)
        } else {
            // Reject the stepsize and propose a smaller one for the current step
-            TryStep::NotYetAccepted(prev_step / (self.factor_c1.min(factor_11 / self.safety_factor)))
+            TryStep::NotYetAccepted(prev_step / self.factor_c1.min(factor_11 / self.safety_factor))
        }
    }
+
+    fn next_step_guess(&self) -> TryStep<T> {
+        self.next_step_guess
+    }
+
+    fn set_next_step_guess(&mut self, guess: TryStep<T>) {
+        self.next_step_guess = guess;
+    }
+
+    fn reset(&mut self) {
+        self.factor_old = 1.0e-4;
+    }
 }

-impl PIController {
+impl<T: Scalar> PIController<T> {
    pub fn new(
        alpha: f64,
        beta: f64,
@@ -83,12 +111,12 @@ impl PIController {
            factor_old: 1.0e-4,
            h_max: h_max.abs(),
            safety_factor,
-            next_step_guess: TryStep::NotYetAccepted(initial_h),
+            next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)),
        }
    }
 }

-impl Default for PIController {
+impl<T: Scalar> Default for PIController<T> {
    fn default() -> Self {
        Self::new(0.17, 0.04, 10.0, 0.2, 100000.0, 0.9, 1e-4)
    }
@@ -111,9 +139,9 @@ impl Default for PIController {
 /// Where:
 /// - ε_n = normalized error estimate at current step
 /// - ε_{n-1} = normalized error estimate at previous step
-/// - β₁ = proportional coefficient (controls reaction to current error)
-/// - β₂ = integral coefficient (controls reaction to error history)
-/// - β₃ = derivative coefficient (controls reaction to error rate of change)
+/// - β₁ (`p_term`) = proportional coefficient (controls reaction to current error)
+/// - β₂ (`i_term`) = integral coefficient (controls reaction to error history)
+/// - β₃ (`d_term`) = derivative coefficient (controls reaction to error rate of change)
 ///
 /// # When to Use
 ///
@@ -138,58 +166,59 @@ impl Default for PIController {
 /// let controller = PIDController::for_order(5);
 /// ```
 #[derive(Debug, Clone, Copy)]
-pub struct PIDController {
+pub struct PIDController<T: Scalar = f64> {
    // PID Coefficients
-    pub beta1: f64,  // Proportional: reaction to current error
-    pub beta2: f64,  // Integral: reaction to error history
-    pub beta3: f64,  // Derivative: reaction to error rate of change
+    pub p_term: f64, // Proportional: reaction to current error
+    pub i_term: f64, // Integral: reaction to error history
+    pub d_term: f64, // Derivative: reaction to error rate of change

    // Constraints
-    pub factor_c1: f64,      // 1/min_factor (maximum step decrease)
-    pub factor_c2: f64,      // 1/max_factor (maximum step increase)
-    pub h_max: f64,          // Maximum allowed step size
-    pub safety_factor: f64,  // Safety factor (typically 0.9)
+    pub factor_c1: f64,     // 1/min_factor (maximum step decrease)
+    pub factor_c2: f64,     // 1/max_factor (maximum step increase)
+    pub h_max: f64,         // Maximum allowed step size
+    pub safety_factor: f64, // Safety factor (typically 0.9)

    // Error history for PID control
-    pub err_old: f64,     // ε_{n-1}: previous step error
-    pub err_older: f64,   // ε_{n-2}: error two steps ago
-    pub h_old: f64,       // h_{n-1}: previous step size
+    pub err_old: f64,   // ε_{n-1}: previous step error
+    pub err_older: f64, // ε_{n-2}: error two steps ago
+    pub h_old: f64,     // h_{n-1}: previous step size

    // Next step guess
-    pub next_step_guess: TryStep,
+    pub next_step_guess: TryStep<T>,
 }

-impl<const D: usize> Controller<D> for PIDController {
+impl<T: Scalar> Controller<T> for PIDController<T> {
    /// Determines if the previously run step was acceptable and computes the next step size
    /// using PID control theory
-    fn determine_step(&mut self, h: f64, err: f64) -> TryStep {
+    fn determine_step(&mut self, h: T, err: f64) -> TryStep<T> {
        // Compute PID control factor
        // For first step or when history isn't available, fall back to simpler control
        let factor = if self.err_old <= 0.0 {
            // First step: use only proportional control
-            let factor_11 = err.powf(self.beta1);
-            self.factor_c2.max(
-                self.factor_c1.min(factor_11 / self.safety_factor)
-            )
+            let factor_11 = err.powf(self.p_term);
+            self.factor_c2
+                .max(self.factor_c1.min(factor_11 / self.safety_factor))
        } else if self.err_older <= 0.0 {
            // Second step: use PI control (proportional + integral)
-            let factor_11 = err.powf(self.beta1);
-            let factor_12 = self.err_old.powf(-self.beta2);
+            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)
-            let factor_11 = err.powf(self.beta1);
-            let factor_12 = self.err_old.powf(-self.beta2);
+            let factor_11 = err.powf(self.p_term);
+            let factor_12 = self.err_old.powf(-self.i_term);
            // Derivative term uses ratio of consecutive step sizes
            let factor_13 = if self.h_old > 0.0 {
-                (h / self.h_old).powf(-self.beta3)
+                (h.re() / self.h_old).powf(-self.d_term)
            } else {
                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),
            )
        };

@@ -199,44 +228,59 @@ impl<const D: usize> Controller<D> for PIDController {

            // 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.h_old = h;
+            self.err_old = err.max(1.0e-4); // Prevent very small values
+            self.h_old = h.re();

            // Apply maximum step size limit
-            if h_next.abs() > self.h_max {
-                h_next = self.h_max.copysign(h_next);
+            if h_next.re().abs() > self.h_max {
+                h_next = T::lift(self.h_max.copysign(h_next.re()));
            }

            TryStep::Accepted(h, h_next)
        } else {
            // Step rejected - propose smaller step
            // Use only proportional control for rejection (more aggressive)
-            let factor_11 = err.powf(self.beta1);
+            let factor_11 = err.powf(self.p_term);
            let h_next = h / (self.factor_c1.min(factor_11 / self.safety_factor));

            // Note: Don't update history on rejection
            TryStep::NotYetAccepted(h_next)
        }
    }
+
+    fn next_step_guess(&self) -> TryStep<T> {
+        self.next_step_guess
+    }
+
+    fn set_next_step_guess(&mut self, guess: TryStep<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 PIDController {
+impl<T: Scalar> PIDController<T> {
    /// Create a new PID controller with custom parameters
    ///
    /// # Arguments
    ///
-    /// * `beta1` - Proportional coefficient (typically 0.3-0.5)
-    /// * `beta2` - Integral coefficient (typically 0.04-0.1)
-    /// * `beta3` - Derivative coefficient (typically 0.01-0.05)
+    /// * `p_term` - Proportional coefficient (typically 0.3-0.5)
+    /// * `i_term` - Integral coefficient (typically 0.04-0.1)
+    /// * `d_term` - Derivative coefficient (typically 0.01-0.05)
    /// * `max_factor` - Maximum step size increase factor (typically 10.0)
    /// * `min_factor` - Maximum step size decrease factor (typically 0.2)
    /// * `h_max` - Maximum allowed step size
    /// * `safety_factor` - Safety factor (typically 0.9)
    /// * `initial_h` - Initial step size guess
+    #[allow(clippy::too_many_arguments)]
    pub fn new(
-        beta1: f64,
-        beta2: f64,
-        beta3: f64,
+        p_term: f64,
+        i_term: f64,
+        d_term: f64,
        max_factor: f64,
        min_factor: f64,
        h_max: f64,
@@ -244,26 +288,26 @@ impl PIDController {
        initial_h: f64,
    ) -> Self {
        Self {
-            beta1,
-            beta2,
-            beta3,
+            p_term,
+            i_term,
+            d_term,
            factor_c1: 1.0 / min_factor,
            factor_c2: 1.0 / max_factor,
            h_max: h_max.abs(),
            safety_factor,
-            err_old: 0.0,      // No history initially
-            err_older: 0.0,    // No history initially
-            h_old: 0.0,        // No history initially
-            next_step_guess: TryStep::NotYetAccepted(initial_h),
+            err_old: 0.0,   // No history initially
+            err_older: 0.0, // No history initially
+            h_old: 0.0,     // No history initially
+            next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)),
        }
    }

    /// Create a PID controller with coefficients scaled for a specific method order
    ///
    /// Uses the Gustafsson coefficients scaled by order:
-    /// - β₁ = 0.49 / (order + 1)
-    /// - β₂ = 0.34 / (order + 1)
-    /// - β₃ = 0.10 / (order + 1)
+    /// - `p_term` = 0.49 / (order + 1)
+    /// - `i_term` = 0.34 / (order + 1)
+    /// - `d_term` = 0.10 / (order + 1)
    ///
    /// # Arguments
    ///
@@ -271,48 +315,39 @@ impl PIDController {
    pub fn for_order(order: usize) -> Self {
        let k_plus_1 = (order + 1) as f64;
        Self::new(
-            0.49 / k_plus_1,  // beta1: proportional
-            0.34 / k_plus_1,  // beta2: integral
-            0.10 / k_plus_1,  // beta3: derivative
-            10.0,              // max_factor
-            0.2,               // min_factor
-            100000.0,          // h_max
-            0.9,               // safety_factor
-            1e-4,              // initial_h
+            0.49 / k_plus_1, // p_term: proportional
+            0.34 / k_plus_1, // i_term: integral
+            0.10 / k_plus_1, // d_term: derivative
+            10.0,            // max_factor
+            0.2,             // min_factor
+            100000.0,        // h_max
+            0.9,             // safety_factor
+            1e-4,            // initial_h
        )
    }
-
-    /// Reset the controller's error history
-    ///
-    /// Useful when switching algorithms or restarting integration
-    pub fn reset(&mut self) {
-        self.err_old = 0.0;
-        self.err_older = 0.0;
-        self.h_old = 0.0;
-    }
 }

-impl Default for PIDController {
+impl<T: Scalar> Default for PIDController<T> {
    /// Default PID controller using conservative coefficients
    ///
    /// Uses conservative PID coefficients that provide stable performance
    /// across a wide range of problems:
-    /// - β₁ = 0.07 (proportional)
-    /// - β₂ = 0.04 (integral)
-    /// - β₃ = 0.01 (derivative)
+    /// - `p_term` = 0.07 (proportional)
+    /// - `i_term` = 0.04 (integral)
+    /// - `d_term` = 0.01 (derivative)
    ///
    /// These values are appropriate for typical ODE methods of order 5-7.
    /// For method-specific tuning, use `PIDController::for_order(order)` instead.
    fn default() -> Self {
        Self::new(
-            0.07,      // beta1 (proportional)
-            0.04,      // beta2 (integral)
-            0.01,      // beta3 (derivative)
-            10.0,      // max_factor
-            0.2,       // min_factor
-            100000.0,  // h_max
-            0.9,       // safety_factor
-            1e-4,      // initial_h
+            0.07,     // p_term (proportional)
+            0.04,     // i_term (integral)
+            0.01,     // d_term (derivative)
+            10.0,     // max_factor
+            0.2,      // min_factor
+            100000.0, // h_max
+            0.9,      // safety_factor
+            1e-4,     // initial_h
        )
    }
 }
@@ -323,7 +358,7 @@ mod tests {

    #[test]
    fn test_pi_controller_creation() {
-        let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 10.0, 0.9, 1e-4);
+        let controller = PIController::<f64>::new(0.17, 0.04, 10.0, 0.2, 10.0, 0.9, 1e-4);

        assert!(controller.alpha == 0.17);
        assert!(controller.beta == 0.04);
@@ -337,11 +372,11 @@ mod tests {

    #[test]
    fn test_pid_controller_creation() {
-        let controller = PIDController::new(0.49, 0.34, 0.10, 10.0, 0.2, 10.0, 0.9, 1e-4);
+        let controller = PIDController::<f64>::new(0.49, 0.34, 0.10, 10.0, 0.2, 10.0, 0.9, 1e-4);

-        assert_eq!(controller.beta1, 0.49);
-        assert_eq!(controller.beta2, 0.34);
-        assert_eq!(controller.beta3, 0.10);
+        assert_eq!(controller.p_term, 0.49);
+        assert_eq!(controller.i_term, 0.34);
+        assert_eq!(controller.d_term, 0.10);
        assert_eq!(controller.factor_c1, 1.0 / 0.2);
        assert_eq!(controller.factor_c2, 1.0 / 10.0);
        assert_eq!(controller.h_max, 10.0);
@@ -353,32 +388,33 @@ mod tests {

    #[test]
    fn test_pid_for_order() {
-        let controller = PIDController::for_order(5);
+        let controller = PIDController::<f64>::for_order(5);

        // For order 5, k+1 = 6
-        assert!((controller.beta1 - 0.49 / 6.0).abs() < 1e-10);
-        assert!((controller.beta2 - 0.34 / 6.0).abs() < 1e-10);
-        assert!((controller.beta3 - 0.10 / 6.0).abs() < 1e-10);
+        assert!((controller.p_term - 0.49 / 6.0).abs() < 1e-10);
+        assert!((controller.i_term - 0.34 / 6.0).abs() < 1e-10);
+        assert!((controller.d_term - 0.10 / 6.0).abs() < 1e-10);
    }

    #[test]
    fn test_pid_default() {
-        let controller = PIDController::default();
+        let controller = PIDController::<f64>::default();

        // Default uses conservative PID coefficients
-        assert_eq!(controller.beta1, 0.07);
-        assert_eq!(controller.beta2, 0.04);
-        assert_eq!(controller.beta3, 0.01);  // True PID with derivative term
+        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
    }

    #[test]
    fn test_pid_reset() {
-        let mut controller = PIDController::default();
-
        // Simulate some history
-        controller.err_old = 0.5;
-        controller.err_older = 0.3;
-        controller.h_old = 0.01;
+        let mut controller = PIDController::<f64> {
+            err_old: 0.5,
+            err_older: 0.3,
+            h_old: 0.01,
+            ..Default::default()
+        };

        controller.reset();

@@ -392,8 +428,8 @@ mod tests {
        // For smooth problems, PI and PID should perform similarly
        // Test with exponential decay: y' = -y

-        let mut pi = PIController::default();
-        let mut pid = PIDController::default();
+        let mut pi = PIController::<f64>::default();
+        let mut pid = PIDController::<f64>::default();

        // Simulate a sequence of small errors (smooth problem)
        let errors = vec![0.8, 0.6, 0.5, 0.45, 0.4, 0.35, 0.3];
@@ -403,8 +439,9 @@ mod tests {
        let mut pid_steps = Vec::new();

        for &err in &errors {
-            let mut pi_result = <PIController as Controller<1>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<1>>::determine_step(&mut pid, h, err);
+            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);

            if pi_result.is_accepted() {
                pi_steps.push(pi_result.extract());
@@ -437,25 +474,25 @@ mod tests {
    #[test]
    fn test_pid_bootstrap() {
        // Test that PID progressively uses P → PI → PID as history builds
-        let mut pid = PIDController::new(0.49, 0.34, 0.10, 10.0, 0.2, 100.0, 0.9, 0.01);
+        let mut pid = PIDController::<f64>::new(0.49, 0.34, 0.10, 10.0, 0.2, 100.0, 0.9, 0.01);

        let h = 0.01;
        let err1 = 0.5;
        let err2 = 0.4;
        let err3 = 0.3;

-        // First step: should use only proportional (beta1)
+        // First step: should use only proportional (p_term)
        assert_eq!(pid.err_old, 0.0);
        assert_eq!(pid.err_older, 0.0);
-        let step1 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err1);
+        let step1 = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err1);
        assert!(step1.is_accepted());

        // After first step, err_old is updated but err_older is still 0
        assert!(pid.err_old > 0.0);
        assert_eq!(pid.err_older, 0.0);

-        // Second step: should use PI (beta1 and beta2)
-        let step2 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err2);
+        // Second step: should use PI (p_term and i_term)
+        let step2 = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err2);
        assert!(step2.is_accepted());

        // After second step, both err_old and err_older should be set
@@ -463,21 +500,21 @@ mod tests {
        assert!(pid.err_older > 0.0);
        assert!(pid.h_old > 0.0);

-        // Third step: should use full PID (beta1, beta2, and beta3)
-        let step3 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err3);
+        // Third step: should use full PID (p_term, i_term, and d_term)
+        let step3 = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err3);
        assert!(step3.is_accepted());
    }

    #[test]
    fn test_pid_step_rejection() {
        // Test that error history is NOT updated after rejection
-        let mut pid = PIDController::default();
+        let mut pid = PIDController::<f64>::default();

        let h = 0.01;

        // First, accept a step to build history
        let err_good = 0.5;
-        let step1 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err_good);
+        let step1 = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err_good);
        assert!(step1.is_accepted());

        let err_old_before = pid.err_old;
@@ -486,7 +523,7 @@ mod tests {

        // Now reject a step with large error
        let err_bad = 2.0;
-        let step2 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err_bad);
+        let step2 = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err_bad);
        assert!(!step2.is_accepted());

        // History should NOT have changed after rejection
@@ -500,9 +537,9 @@ mod tests {
        // Test on oscillatory error pattern (simulating difficult problem)
        // True PID (with derivative term) should provide smoother step size evolution

-        let mut pi = PIController::default();
-        // Use actual PID with non-zero beta3 (Gustafsson coefficients for order 5)
-        let mut pid = PIDController::for_order(5);
+        let mut pi = PIController::<f64>::default();
+        // Use actual PID with non-zero d_term (Gustafsson coefficients for order 5)
+        let mut pid = PIDController::<f64>::for_order(5);

        // Simulate oscillatory error pattern
        let errors = vec![0.8, 0.3, 0.9, 0.2, 0.85, 0.25, 0.8, 0.3];
@@ -515,8 +552,9 @@ mod tests {
        let mut pid_h_prev = h;

        for &err in &errors {
-            let mut pi_result = <PIController as Controller<1>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<1>>::determine_step(&mut pid, h, err);
+            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);

            if pi_result.is_accepted() {
                let pi_h_next = pi_result.reset().unwrap().extract();
@@ -535,11 +573,8 @@ mod tests {

        // Compute standard deviation of log step size ratios
        let pi_mean: f64 = pi_ratios.iter().sum::<f64>() / pi_ratios.len() as f64;
-        let pi_variance: f64 = pi_ratios
-            .iter()
-            .map(|r| (r - pi_mean).powi(2))
-            .sum::<f64>()
-            / pi_ratios.len() as f64;
+        let pi_variance: f64 =
+            pi_ratios.iter().map(|r| (r - pi_mean).powi(2)).sum::<f64>() / pi_ratios.len() as f64;
        let pi_std = pi_variance.sqrt();

        let pid_mean: f64 = pid_ratios.iter().sum::<f64>() / pid_ratios.len() as f64;
@@ -550,7 +585,7 @@ mod tests {
            / pid_ratios.len() as f64;
        let pid_std = pid_variance.sqrt();

-        // With true PID (non-zero beta3), we expect similar or better stability
+        // With true PID (non-zero d_term), we expect similar or better stability
        // Allow some tolerance since this is a simple synthetic test
        assert!(
            pid_std <= pi_std * 1.1,
@@ -1,7 +1,11 @@
 use nalgebra::SVector;

 use super::super::ode::ODE;
-use super::Integrator;
+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> {
@@ -57,6 +61,15 @@ pub struct BS3<const D: usize> {
    r_tol: f64,
 }

+impl<const D: usize> Default for BS3<D>
+where
+    BS3<D>: Integrator<D>,
+{
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
 impl<const D: usize> BS3<D>
 where
    BS3<D>: Integrator<D>,
@@ -70,21 +83,31 @@ where
            r_tol: 1e-8,
        }
    }
+}

-    /// Set absolute tolerance (same value for all components)
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for BS3<D>
+where
+    BS3<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
        self.a_tol = SVector::<f64, D>::from_element(a_tol);
        self
    }

-    /// Set absolute tolerance (different value per component)
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
        self.a_tol = a_tol;
        self
    }

-    /// Set relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
+    fn r_tol(mut self, r_tol: f64) -> Self {
        self.r_tol = r_tol;
        self
    }
@@ -97,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]
@@ -107,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]
@@ -115,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]
@@ -123,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]
@@ -135,64 +168,95 @@ impl<'a, const D: usize> Integrator<D> for BS3<D>
 where
    BS3<D>: BS3Integrator<'a>,
 {
-    const ORDER: usize = 3;
-    const STAGES: usize = 4;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        3
+    }

-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
+    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<f64, D>> = vec![SVector::<f64, 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);

        // Stage 2: k2 = f(t + c[1]*h, y + h*a[1,0]*k1)
-        let y2 = ode.y + h * Self::A[0] * k[0];
-        k[1] = (ode.f)(ode.t + Self::C[1] * h, y2, &ode.params);
+        let y2 = ode.y + k[0] * (h * Self::A[0]);
+        k[1] = (ode.f)(ode.t + h * Self::C[1], y2, &ode.params);

        // Stage 3: k3 = f(t + c[2]*h, y + h*(a[2,0]*k1 + a[2,1]*k2))
-        let y3 = ode.y + h * (Self::A[1] * k[0] + Self::A[2] * k[1]);
-        k[2] = (ode.f)(ode.t + Self::C[2] * h, y3, &ode.params);
+        let y3 = ode.y + (k[0] * T::lift(Self::A[1]) + k[1] * T::lift(Self::A[2])) * h;
+        k[2] = (ode.f)(ode.t + h * Self::C[2], y3, &ode.params);

        // 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 + h * (Self::A[3] * k[0] + Self::A[4] * k[1] + Self::A[5] * k[2]);
-        k[3] = (ode.f)(ode.t + Self::C[3] * h, y4, &ode.params);
+        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]))
+                * h;
+        k[3] = (ode.f)(ode.t + h * Self::C[3], y4, &ode.params);

        // Compute 3rd order solution
-        let next_y = ode.y + h * (Self::B[0] * k[0] + Self::B[1] * k[1] + Self::B[2] * k[2] + Self::B[3] * k[3]);
+        let next_y = ode.y
+            + (k[0] * T::lift(Self::B[0])
+                + k[1] * T::lift(Self::B[1])
+                + k[2] * T::lift(Self::B[2])
+                + k[3] * T::lift(Self::B[3]))
+                * h;

        // Compute error estimate (difference between 3rd and 2nd order solutions)
-        let err = h * (Self::B_ERROR[0] * k[0] + Self::B_ERROR[1] * k[1] + Self::B_ERROR[2] * k[2] + Self::B_ERROR[3] * k[3]);
+        let err = (k[0] * T::lift(Self::B_ERROR[0])
+            + k[1] * T::lift(Self::B_ERROR[1])
+            + k[2] * T::lift(Self::B_ERROR[2])
+            + k[3] * T::lift(Self::B_ERROR[3]))
+            * h;

-        // Compute error norm scaled by tolerance
-        let tol = self.a_tol + ode.y.abs() * self.r_tol;
-        let error_norm = (err.component_div(&tol)).norm();
+        // Compute error norm scaled by tolerance. Step control is a real-valued
+        // decision, so each component drops to the reals as it is accumulated —
+        // one pass, rather than building real copies of `err` and `y` first.
+        let mut error_norm_sq = 0.0;
+        for i in 0..D {
+            let tol = self.a_tol[i] + ode.y[i].re().abs() * self.r_tol;
+            let scaled = err[i].re() / tol;
+            error_norm_sq += scaled * scaled;
+        }
+        let error_norm = error_norm_sq.sqrt();

        // Store coefficients for dense output (cubic Hermite interpolation)
        // BS3 uses standard cubic Hermite interpolation with derivatives at endpoints
        // Store: y0, y1, f0=k[0], f1=k[3] (FSAL)
        let dense_coeffs = vec![
-            ode.y,           // y0 at start of step
-            next_y,          // y1 at end of step
-            k[0],            // f(t0, y0) - derivative at start
-            k[3],            // f(t1, y1) - derivative at end (FSAL)
+            ode.y,  // y0 at start of step
+            next_y, // y1 at end of step
+            k[0],   // f(t0, y0) - derivative at start
+            k[3],   // f(t1, y1) - derivative at end (FSAL)
        ];

-        (next_y, Some(error_norm), Some(dense_coeffs))
+        StepResult {
+            y: next_y,
+            err: Some(error_norm),
+            dense: Some(dense_coeffs),
+        }
    }

-    fn interpolate(
+    fn interpolate<T: Scalar>(
        &self,
-        t_start: f64,
-        t_end: f64,
-        dense: &[SVector<f64, D>],
-        t: f64,
-    ) -> SVector<f64, D> {
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D> {
        // Compute interpolation parameter θ ∈ [0, 1]
        let theta = (t - t_start) / (t_end - t_start);
        let h = t_end - t_start;
@@ -215,14 +279,14 @@ where
        let f1 = &dense[3];

        let theta2 = theta * theta;
-        let one_minus_theta = 1.0 - theta;
+        let one_minus_theta = T::lift(1.0) - theta;
        let one_minus_theta2 = one_minus_theta * one_minus_theta;

        // Apply cubic Hermite interpolation formula
-        (1.0 + 2.0 * theta) * one_minus_theta2 * y0
-            + theta2 * (3.0 - 2.0 * theta) * y1
-            + theta * one_minus_theta2 * h * f0
-            + theta2 * (theta - 1.0) * h * f1
+        y0 * ((theta * 2.0 + 1.0) * one_minus_theta2)
+            + y1 * (theta2 * (-(theta * 2.0 - 3.0)))
+            + f0 * (theta * one_minus_theta2 * h)
+            + f1 * (theta2 * (theta - 1.0) * h)
    }
 }

@@ -234,11 +298,11 @@ mod tests {

    #[test]
    fn test_bs3_creation() {
-        let _bs3: BS3<1> = BS3::new();
-        assert_eq!(BS3::<1>::ORDER, 3);
-        assert_eq!(BS3::<1>::STAGES, 4);
-        assert!(BS3::<1>::ADAPTIVE);
-        assert!(BS3::<1>::DENSE);
+        let bs3: BS3<1> = BS3::new();
+        assert_eq!(bs3.order(), 3);
+        assert_eq!(bs3.stages(), 4);
+        assert!(bs3.adaptive());
+        assert!(bs3.dense());
    }

    #[test]
@@ -252,9 +316,10 @@ 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 (y_next, err, dense) = bs3.step(&ode, h);
+        let __step = bs3.step(&ode, h);
+        let (y_next, err, dense) = (__step.y, __step.err, __step.dense);

        // At t=0.001, exact solution is e^0.001 ≈ 1.0010005001667084
        let exact = (0.001_f64).exp();
@@ -281,9 +346,10 @@ 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 (_y_next, _err, dense) = bs3.step(&ode, h);
+        let __step = bs3.step(&ode, h);
+        let (_y_next, _err, dense) = (__step.y, __step.err, __step.dense);
        let dense = dense.unwrap();

        // Interpolate at midpoint
@@ -312,7 +378,8 @@ mod tests {
        // Take 100 steps to reach t = 1.0
        let mut ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
        for _ in 0..100 {
-            let (y_new, _, _) = bs3.step(&ode, h);
+            let __step = bs3.step(&ode, h);
+            let (y_new, _, _) = (__step.y, __step.err, __step.dense);
            ode.y = y_new;
            ode.t += h;
        }
@@ -345,7 +412,8 @@ mod tests {

            // Take steps until we reach t_end
            while ode.t < t_end - 1e-10 {
-                let (y_new, _, _) = bs3.step(&ode, h);
+                let __step = bs3.step(&ode, h);
+                let (y_new, _, _) = (__step.y, __step.err, __step.dense);
                ode.y = y_new;
                ode.t += h;
            }
@@ -389,7 +457,8 @@ mod tests {

        // First step
        let ode1 = ODE::new(&derivative, 0.0, 1.0, y0, ());
-        let (y_new1, _, dense1) = bs3.step(&ode1, h);
+        let __step = bs3.step(&ode1, h);
+        let (y_new1, _, dense1) = (__step.y, __step.err, __step.dense);
        let dense1 = dense1.unwrap();

        // Extract f1 from first step (derivative at end of step)
@@ -397,7 +466,8 @@ mod tests {

        // Second step starts where first ended
        let ode2 = ODE::new(&derivative, h, 1.0, y_new1, ());
-        let (_, _, dense2) = bs3.step(&ode2, h);
+        let __step = bs3.step(&ode2, h);
+        let (_, _, dense2) = (__step.y, __step.err, __step.dense);
        let dense2 = dense2.unwrap();

        // Extract f0 from second step (derivative at start of step)
@@ -1,7 +1,12 @@
 use nalgebra::SVector;

-use super::super::ode::ODE;
-use super::Integrator;
+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> {
@@ -11,9 +16,54 @@ pub trait DormandPrinceIntegrator<'a> {
    const D: &'a [f64];
 }

+/// Dormand-Prince 5(4) method
+///
+/// A 5th order explicit Runge-Kutta method with an embedded 4th order method for
+/// error estimation. This is the general-purpose workhorse of the crate, and the
+/// same method as MATLAB's `ode45` and `DifferentialEquations.jl`'s `DP5`.
+///
+/// # Characteristics
+/// - Order: 5(4) - 5th order solution with 4th order error estimate
+/// - Stages: 7
+/// - FSAL: Yes (First Same As Last - the 7th stage is the next step's first,
+///   though this implementation does not yet reuse it)
+/// - Adaptive: Yes
+/// - Dense output: 4th order interpolation
+///
+/// # When to use DormandPrince45
+/// - The default choice when nothing about the problem suggests otherwise
+/// - Problems requiring general accuracy (rtol ~ 1e-6 to 1e-8)
+/// - Non-stiff problems
+///
+/// At tighter tolerances than about 1e-8, prefer [`Vern7`](super::vern7::Vern7):
+/// a higher order method reaches the same accuracy in far fewer steps.
+///
+/// # Example
+/// ```rust
+/// use ordinary_diffeq::prelude::*;
+/// use nalgebra::Vector1;
+///
+/// let params = ();
+/// fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+///     Vector1::new(-y[0])
+/// }
+///
+/// let y0 = Vector1::new(1.0);
+/// let ode = ODE::new(&derivative, 0.0, 5.0, y0, ());
+/// let dp45 = DormandPrince45::new().a_tol(1e-8).r_tol(1e-8);
+/// let controller = PIController::default();
+///
+/// let mut problem = Problem::new(ode, dp45, controller);
+/// let solution = problem.solve();
+/// ```
+///
+/// # References
+/// - Dormand, J.R. and Prince, P.J. (1980), "A family of embedded Runge-Kutta
+///   formulae", Journal of Computational and Applied Mathematics, Vol. 6, No. 1,
+///   pp. 19-26
 #[derive(Debug, Clone, Copy)]
 pub struct DormandPrince45<const D: usize> {
-    a_tol: SVector<f64,D>,
+    a_tol: SVector<f64, D>,
    r_tol: f64,
 }

@@ -22,16 +72,53 @@ where
    DormandPrince45<D>: Integrator<D>,
 {
    pub fn new() -> Self {
-        Self { a_tol: SVector::<f64,D>::from_element(1e-8), r_tol: 1e-8 }
+        Self {
+            a_tol: SVector::<f64, D>::from_element(1e-8),
+            r_tol: 1e-8,
+        }
    }
-    pub fn a_tol(&mut self, a_tol: f64) -> Self {
-        Self { a_tol: SVector::<f64,D>::from_element(a_tol), r_tol: self.r_tol }
+}
+
+impl<const D: usize> Tolerances<D> for DormandPrince45<D>
+where
+    DormandPrince45<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
    }
-    pub fn a_tol_full(&mut self, a_tol: SVector::<f64,D>) -> Self {
-        Self { a_tol, r_tol: self.r_tol }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
    }
-    pub fn r_tol(&mut self, r_tol: f64) -> Self {
-        Self { a_tol: self.a_tol, 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,
+        }
+    }
+
+    fn a_tol_full(self, a_tol: SVector<f64, D>) -> Self {
+        Self {
+            a_tol,
+            r_tol: self.r_tol,
+        }
+    }
+
+    fn r_tol(self, r_tol: f64) -> Self {
+        Self {
+            a_tol: self.a_tol,
+            r_tol,
+        }
+    }
+}
+
+impl<const D: usize> Default for DormandPrince45<D>
+where
+    DormandPrince45<D>: Integrator<D>,
+{
+    fn default() -> Self {
+        Self::new()
    }
 }

@@ -91,56 +178,78 @@ impl<'a, const D: usize> Integrator<D> for DormandPrince45<D>
 where
    DormandPrince45<D>: DormandPrinceIntegrator<'a>,
 {
-    const ORDER: usize = 5;
-    const STAGES: usize = 7;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        5
+    }

-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
-        let mut k: Vec<SVector<f64, D>> = vec![SVector::<f64, D>::zeros(); Self::STAGES];
+    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(); STAGES];
        let mut next_y = ode.y;
-        let mut err = SVector::<f64, D>::zeros();
-        let mut rcont5 = SVector::<f64, D>::zeros();
+        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] * Self::B[0] * h;
-        err += k[0] * (Self::B[0] - Self::B[Self::STAGES]) * h;
+        next_y += k[0] * T::lift(Self::B[0]) * h;
+        err += k[0] * T::lift(Self::B[0] - Self::B[STAGES]) * h;
        let rcont1 = ode.y;
-        rcont5 += k[0] * h * Self::D[0];
+        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::<f64, D>::zeros();
+            let mut y_term = SVector::<T, D>::zeros();
            for (j, item) in k.iter().enumerate().take(i) {
-                y_term += item * Self::A[(i * (i - 1)) / 2 + j];
+                y_term += item * T::lift(Self::A[(i * (i - 1)) / 2 + j]);
            }
-            k[i] = (ode.f)(ode.t + 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 * Self::B[i];
-            err += k[i] * (Self::B[i] - Self::B[i + Self::STAGES]) * h;
-            rcont5 += k[i] * h * Self::D[i];
+            next_y += k[i] * h * T::lift(Self::B[i]);
+            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 = h * k[0] - rcont2;
-        let rcont4 = rcont2 - k[Self::STAGES - 1] * h - rcont3;
-        let tol = self.a_tol + ode.y * self.r_tol;
+        let rcont3 = k[0] * h - rcont2;
+        let rcont4 = rcont2 - k[STAGES - 1] * h - rcont3;
        let rcont = vec![rcont1, rcont2, rcont3, rcont4, rcont5];
-        (next_y, Some((err.component_div(&tol)).norm()), Some(rcont))
+        // Just take the real component of the error
+        let mut error_norm_sq = 0.0;
+        for i in 0..D {
+            let tol = self.a_tol[i] + ode.y[i].re() * self.r_tol;
+            let scaled = err[i].re() / tol;
+            error_norm_sq += scaled * scaled;
+        }
+        StepResult {
+            y: next_y,
+            err: Some(error_norm_sq.sqrt()),
+            dense: Some(rcont),
+        }
    }
-    fn interpolate(
+
+    fn interpolate<T: Scalar>(
        &self,
-        t_start: f64,
-        t_end: f64,
-        dense: &[SVector<f64, D>],
-        t: f64,
-    ) -> SVector<f64, D> {
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D> {
        let s = (t - t_start) / (t_end - t_start);
-        let s1 = 1.0 - s;
+        let s1 = T::one() - s;
        dense[0] + (dense[1] + (dense[2] + (dense[3] + dense[4] * s1) * s) * s1) * s
    }
 }
@@ -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));
+    }
+}
@@ -1,43 +1,84 @@
 use nalgebra::SVector;

-use super::ode::ODE;
+use crate::ode::ODE;
+use crate::scalar::Scalar;

 pub mod bs3;
 pub mod dormand_prince;
+pub mod method;
 pub mod rosenbrock;
 pub mod vern7;

-/// Integrator Trait
-pub trait Integrator<const D: usize> {
-    const ORDER: usize;
-    const STAGES: usize;
-    const ADAPTIVE: bool;
-    const DENSE: bool;
+/// Result of a `Integrator::step`
+#[derive(Debug, Clone, PartialEq)]
+pub struct StepResult<T: Scalar, const D: usize> {
+    /// The state at the end of the step.
+    pub y: SVector<T, D>,

-    /// Number of main stages stored in dense output (default: same as STAGES)
-    const MAIN_STAGES: usize = Self::STAGES;
+    /// The local error, normalized by the integrator's tolerances, or `None`
+    /// from a fixed-step method. A value at or below `1.0` means the step is
+    /// acceptable.
+    pub err: Option<f64>,
+
+    /// Coefficients for interpolating inside the step, to be handed back to
+    /// [`Integrator::interpolate`], or `None` from a method without dense
+    /// output.
+    ///
+    /// For a method with lazy dense output (`extra_stages() > 0`) the layout is
+    /// `[y_start, y_end, k_1, ..., k_main_stages]`, and the extra stages are
+    /// 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> {
+    /// The order of the method's solution.
+    fn order(&self) -> usize;
+
+    /// The number of stages the method evaluates per step.
+    fn stages(&self) -> usize;
+
+    /// Whether the method estimates its local error, and so can be stepped
+    /// adaptively.
+    fn adaptive(&self) -> bool;
+
+    /// Whether the method produces dense output.
+    fn dense(&self) -> bool;
+
+    /// Number of main stages stored in dense output (default: same as `stages`)
+    fn main_stages(&self) -> usize {
+        self.stages()
+    }

    /// Number of extra stages for full-order dense output (default: 0, no extra stages)
-    const EXTRA_STAGES: usize = 0;
+    fn extra_stages(&self) -> usize {
+        0
+    }

    /// Total stages when full dense output is computed
-    const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES;
+    fn total_dense_stages(&self) -> usize {
+        self.main_stages() + self.extra_stages()
+    }

-    /// Returns a new y value, then possibly an error value, and possibly a dense output
-    /// coefficient set
-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>);
+    /// 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>;

-    fn interpolate(
+    /// Evaluates the solution at `t` inside a step that ran from `t_start` to
+    /// `t_end`, using the `dense` coefficients that step produced.
+    fn interpolate<T: Scalar>(
        &self,
-        t_start: f64,
-        t_end: f64,
-        dense: &[SVector<f64, D>],
-        t: f64,
-    ) -> SVector<f64, D>;
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D>;

    /// Compute extra stages for full-order dense output (lazy computation).
    ///
@@ -56,19 +97,59 @@ pub trait Integrator<const D: usize> {
    /// # Returns
    ///
    /// Vector of extra k-stages (empty for most integrators)
-    fn compute_extra_stages<P>(
+    fn compute_extra_stages<T: Scalar, P>(
        &self,
-        _ode: &ODE<D, P>,
-        _t_start: f64,
-        _y_start: SVector<f64, D>,
-        _h: f64,
-        _main_stages: &[SVector<f64, D>],
-    ) -> Vec<SVector<f64, D>> {
+        _ode: &ODE<T, D, P>,
+        _t_start: T,
+        _y_start: SVector<T, D>,
+        _h: T,
+        _main_stages: &[SVector<T, D>],
+    ) -> Vec<SVector<T, D>> {
        // Default implementation: no extra stages needed
        Vec::new()
    }
 }

+/// 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::*;
@@ -79,7 +160,7 @@ mod tests {
    #[test]
    fn test_dopri() {
        type Params = ();
-        fn derivative(_t: f64, y: Vector3<f64>, _p: &Params) -> Vector3<f64> {
+        fn derivative<T: Scalar>(_t: T, y: Vector3<T>, _p: &Params) -> Vector3<T> {
            y
        }

@@ -92,11 +173,109 @@ mod tests {
        // and also that the error seems reasonable
        let step = 0.001;
        while ode.t < ode.t_end {
-            let (new_y, err, _) = dp45.step(&ode, step);
-            ode.y = new_y;
+            let result = dp45.step(&ode, step);
+            ode.y = result.y;
            ode.t += step;
            assert_relative_eq!(ode.y[0], ode.t.exp(), max_relative = 0.01);
-            assert!(err.unwrap() < 1.0);
+            assert!(result.err.unwrap() < 1.0);
+        }
+    }
+}
+
+#[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"
+            );
        }
    }
 }
@@ -1,7 +1,8 @@
 use nalgebra::{SMatrix, SVector};

-use super::super::ode::ODE;
-use super::Integrator;
+use super::super::ode::{ProblemFunction, ODE};
+use super::{Integrator, StepResult, Tolerances};
+use crate::scalar::Scalar;

 /// Strategy for when to update the Jacobian matrix
 #[derive(Debug, Clone, Copy)]
@@ -27,28 +28,31 @@ impl Default for JacobianUpdateStrategy {

 /// Compute the Jacobian matrix ∂f/∂y using forward finite differences
 ///
-/// For a system y' = f(t, y), computes the D×D Jacobian matrix J where J[i,j] = ∂f_i/∂y_j
+/// For a system `y' = f(t, y)`, computes the D×D Jacobian matrix `J` where
+/// `J[i,j] = ∂f_i/∂y_j`
 ///
-/// Uses forward differences: J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε
-/// where ε = √(machine_epsilon) * max(|y[j]|, 1.0)
-pub fn compute_jacobian<const D: usize, P>(
-    f: &dyn Fn(f64, SVector<f64, D>, &P) -> SVector<f64, D>,
-    t: f64,
-    y: SVector<f64, D>,
+/// Uses forward differences: `J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε`
+/// where `ε = √(machine_epsilon) * max(|y[j]|, 1.0)`
+pub fn compute_jacobian<T: Scalar, const D: usize, P>(
+    f: ProblemFunction<'_, T, D, P>,
+    t: T,
+    y: SVector<T, D>,
    params: &P,
-) -> SMatrix<f64, D, D> {
+) -> SMatrix<T, D, D> {
    let sqrt_eps = f64::EPSILON.sqrt();
    let f_y = f(t, y, params);
-    let mut jacobian = SMatrix::<f64, D, D>::zeros();
+    let mut jacobian = SMatrix::<T, D, D>::zeros();

    // Compute each column of the Jacobian by perturbing y[j]
    for j in 0..D {
-        // Choose epsilon based on the magnitude of y[j]
-        let epsilon = sqrt_eps * y[j].abs().max(1.0);
+        // Choose epsilon based on the magnitude of y[j]. The perturbation is a
+        // numerical-method parameter, not part of the solution, so it is chosen
+        // in the reals.
+        let epsilon = sqrt_eps * y[j].re().abs().max(1.0);

        // Perturb y in the j-th direction
        let mut y_perturbed = y;
-        y_perturbed[j] += epsilon;
+        y_perturbed[j] += T::lift(epsilon);

        // Evaluate f at perturbed point
        let f_perturbed = f(t, y_perturbed, params);
@@ -131,8 +135,8 @@ pub fn compute_jacobian<const D: usize, P>(
 /// let rosenbrock = Rosenbrock23::new();
 ///
 /// // Take a single step
-/// let (y_next, error, _dense) = rosenbrock.step(&ode, 0.1);
-/// assert!((y_next[0] - 0.905).abs() < 0.01);
+/// let step = rosenbrock.step(&ode, 0.1);
+/// assert!((step.y[0] - 0.905).abs() < 0.01);
 /// ```
 #[derive(Debug, Clone, Copy)]
 pub struct Rosenbrock23<const D: usize> {
@@ -140,8 +144,8 @@ pub struct Rosenbrock23<const D: usize> {
    c32: f64,
    /// Coefficient d = 1/(2+√2) ≈ 0.29289321881345254 (this is gamma!)
    d: f64,
-    /// Absolute tolerance for error estimation
-    a_tol: f64,
+    /// Absolute tolerance for error estimation, per component
+    a_tol: SVector<f64, D>,
    /// Relative tolerance for error estimation
    r_tol: f64,
    /// Strategy for updating the Jacobian
@@ -162,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,
@@ -172,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;
@@ -206,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()
@@ -213,86 +230,76 @@ impl<const D: usize> Default for Rosenbrock23<D> {
 }

 impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
-    const ORDER: usize = 2;
-    const STAGES: usize = 2;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        2
+    }

-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
+    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;
        let uprev = ode.y;

        // Compute Jacobian
        let j = compute_jacobian(&ode.f, t, uprev, &ode.params);

-        // Julia: dtγ = dt * d
        let dtgamma = h * self.d;

        // Form W = I - dtγ*J
-        let w = SMatrix::<f64, D, D>::identity() - dtgamma * j;
+        let w = SMatrix::<T, D, D>::identity() - j * dtgamma;
        let w_inv = w.try_inverse().expect("W matrix is singular");

        // Evaluate fsalfirst = f(uprev)
        let fsalfirst = (ode.f)(t, uprev, &ode.params);

        // Stage 1: Solve W * k₁ = f(y) where k₁ is a derivative estimate
-        // Julia stores derivatives in k, not displacements
        let k1 = w_inv * fsalfirst;

        // Stage 2: u = uprev + dt/2 * k₁
-        // Julia line 69
        let dto2 = h / 2.0;
-        let u = uprev + dto2 * k1;
+        let u = uprev + k1 * dto2;

        // Evaluate f₁ = f(u, t + dt/2)
-        // Julia line 71
        let f1 = (ode.f)(t + dto2, u, &ode.params);

        // Stage 2: W * k₂ = f₁ - k₁ + J*k₁
-        // Julia line 80: linsolve_tmp = f₁ - tmp (where tmp = k₁)
-        // This is equivalent to: W * k₂ = f₁ - k₁
-        //   => (I - dtγ*J) * k₂ = f₁ - k₁
-        //   => k₂ = (I - dtγ*J)^{-1} * (f₁ - k₁)
-        // But actually, maybe the RHS should be scaled differently. Let me try: W * k₂ = f₁ + J*k₁
-        // Since W = I - dtγ*J, we have W*k₂ - I*k₂ = -dtγ*J*k₂, so if RHS = f₁ + J*k₁...
-        // Actually, let's just implement exactly what Julia does:
        let rhs2 = f1 - k1;
        let k2_temp = w_inv * rhs2;
-        // Julia then does: k₂ = k₂_temp * neginvdtγ + k₁
-        // But neginvdtγ = -1/(dtγ), which would give huge values.
-        // Let me try without that scaling:
        let k2 = k2_temp + k1;

        // Solution: u = uprev + dt * k₂
-        // Julia line 89
-        let u_final = uprev + h * k2;
+        let u_final = uprev + k2 * h;

        // Error estimation
        // Evaluate fsallast = f(u_final, t + dt)
-        // Julia line 94
        let fsallast = (ode.f)(t + h, u_final, &ode.params);

-        // Julia line 98-99: linsolve_tmp = fsallast - c₃₂*(k₂ - f₁) - 2*(k₁ - fsalfirst) + dt*dT
        // Ignoring dT (time derivative) for autonomous systems
-        let linsolve_tmp3 = fsallast - self.c32 * (k2 - f1) - 2.0 * (k1 - fsalfirst);
+        let linsolve_tmp3 =
+            fsallast - (k2 - f1) * T::lift(self.c32) - (k1 - fsalfirst) * T::lift(2.0);

        // Stage 3 for error estimation: W * k₃ = linsolve_tmp3
        let k3 = w_inv * linsolve_tmp3;

        // Error: dt/6 * (k₁ - 2*k₂ + k₃)
-        // Julia line 115
        let dto6 = h / 6.0;
-        let error_vec = dto6 * (k1 - 2.0 * k2 + k3);
+        let error_vec = (k1 - k2 * T::lift(2.0) + k3) * dto6;

        // Compute scalar error estimate using weighted norm
        let mut error_sum = 0.0;
        for i in 0..D {
-            let scale = self.a_tol + self.r_tol * uprev[i].abs().max(u_final[i].abs());
-            let weighted_error = error_vec[i] / scale;
+            let scale = self.a_tol[i] + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs());
+            let weighted_error = error_vec[i].re() / scale;
            error_sum += weighted_error * weighted_error;
        }
        let error = (error_sum / D as f64).sqrt();
@@ -300,16 +307,20 @@ impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
        // Dense output: store k₁ and k₂
        let dense = Some(vec![k1, k2]);

-        (u_final, Some(error), dense)
+        StepResult {
+            y: u_final,
+            err: Some(error),
+            dense,
+        }
    }

-    fn interpolate(
+    fn interpolate<T: Scalar>(
        &self,
-        t_start: f64,
-        t_end: f64,
-        dense: &[SVector<f64, D>],
-        t: f64,
-    ) -> SVector<f64, D> {
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D> {
        // Second order interpolation using k₁ and k₂
        // For Rosenbrock methods, we use a simple Hermite interpolation
        let k1 = dense[0];
@@ -323,7 +334,7 @@ impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
        // y(t) ≈ y_n + θ*h*((1-θ)*k₁ + θ*k₂)
        // But we don't have y_n stored, so we just return the stage interpolation
        // This is a simple linear interpolation of the derivative
-        (1.0 - theta) * k1 + theta * k2
+        k1 * (T::lift(1.0) - theta) + k2 * theta
    }
 }

@@ -383,7 +394,8 @@ mod tests {
        let ode = ODE::new(&derivative, 0.0, 0.1, y0, ());
        let rb23 = Rosenbrock23::new();

-        let (y_next, err, _) = rb23.step(&ode, 0.1);
+        let __step = rb23.step(&ode, 0.1);
+        let (y_next, err, _) = (__step.y, __step.err, __step.dense);

        let analytical = (-0.1_f64).exp();
        println!("Computed: {}, Analytical: {}", y_next[0], analytical);
@@ -416,7 +428,8 @@ mod tests {
            let rb23 = Rosenbrock23::new();

            while ode.t < t_end - 1e-10 {
-                let (y_next, _, _) = rb23.step(&ode, h);
+                let __step = rb23.step(&ode, h);
+                let (y_next, _, _) = (__step.y, __step.err, __step.dense);
                ode.y = y_next;
                ode.t += h;
            }
@@ -431,12 +444,21 @@ mod tests {
        // So log(error) = 2*log(h) + const
        // Slope should be approximately 2
        for i in 0..errors.len() - 1 {
-            let rate =
-                (errors[i].log10() - errors[i + 1].log10()) / (step_sizes[i].log10() - step_sizes[i + 1].log10());
-            println!("Convergence rate between h={} and h={}: {}", step_sizes[i], step_sizes[i+1], rate);
+            let rate = (errors[i].log10() - errors[i + 1].log10())
+                / (step_sizes[i].log10() - step_sizes[i + 1].log10());
+            println!(
+                "Convergence rate between h={} and h={}: {}",
+                step_sizes[i],
+                step_sizes[i + 1],
+                rate
+            );

            // Should be close to 2 for a 2nd order method
-            assert!(rate > 1.8 && rate < 2.2, "Convergence rate {} is not close to 2", rate);
+            assert!(
+                rate > 1.8 && rate < 2.2,
+                "Convergence rate {} is not close to 2",
+                rate
+            );
        }
    }

@@ -454,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 {
@@ -465,7 +487,8 @@ mod tests {
            let rb23 = Rosenbrock23::new();

            for _ in 0..n_steps {
-                let (y_next, _, _) = rb23.step(&ode, actual_h);
+                let __step = rb23.step(&ode, actual_h);
+                let (y_next, _, _) = (__step.y, __step.err, __step.dense);
                ode.y = y_next;
                ode.t += actual_h;
            }
@@ -491,13 +514,15 @@ mod tests {
        println!("Orders per step refinement: {:?}", orders);

        // Julia tests: @test sim.𝒪est[:final]≈2 atol=0.2
-        assert!((avg_order - 2.0).abs() < 0.2,
-                "Convergence order {:.2} not within 0.2 of expected order 2", avg_order);
+        assert!(
+            (avg_order - 2.0).abs() < 0.2,
+            "Convergence order {:.2} not within 0.2 of expected order 2",
+            avg_order
+        );
    }

    #[test]
    fn test_rosenbrock23_adaptive_solve() {
-        // Julia test: sol = solve(prob, Rosenbrock23()); @test length(sol) < 20
        // This tests that the adaptive solver can efficiently solve prob_ode_linear
        use crate::controller::PIController;
        use crate::problem::Problem;
@@ -516,12 +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()
+        );

-        // Julia test: @test length(sol) < 20
-        assert!(solution.states.len() < 20,
-                "Adaptive solve should complete in less than 20 steps, got {}",
-                solution.states.len());
+        assert!(
+            solution.states.len() < 20,
+            "Adaptive solve should complete in less than 20 steps, got {}",
+            solution.states.len()
+        );

        // Verify final value is accurate
        let analytical = 0.5 * (1.01_f64 * 1.0).exp();
@@ -1,7 +1,15 @@
 use nalgebra::SVector;

 use super::super::ode::ODE;
-use super::Integrator;
+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> {
@@ -69,6 +77,15 @@ pub struct Vern7<const D: usize> {
    r_tol: f64,
 }

+impl<const D: usize> Default for Vern7<D>
+where
+    Vern7<D>: Integrator<D>,
+{
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
 impl<const D: usize> Vern7<D>
 where
    Vern7<D>: Integrator<D>,
@@ -82,21 +99,31 @@ where
            r_tol: 1e-8,
        }
    }
+}

-    /// Set absolute tolerance (same value for all components)
-    pub fn a_tol(mut self, a_tol: f64) -> Self {
+impl<const D: usize> Tolerances<D> for Vern7<D>
+where
+    Vern7<D>: Integrator<D>,
+{
+    fn absolute_tolerance(&self) -> SVector<f64, D> {
+        self.a_tol
+    }
+
+    fn relative_tolerance(&self) -> f64 {
+        self.r_tol
+    }
+
+    fn a_tol(mut self, a_tol: f64) -> Self {
        self.a_tol = SVector::<f64, D>::from_element(a_tol);
        self
    }

-    /// Set absolute tolerance (different value per component)
-    pub fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
+    fn a_tol_full(mut self, a_tol: SVector<f64, D>) -> Self {
        self.a_tol = a_tol;
        self
    }

-    /// Set relative tolerance
-    pub fn r_tol(mut self, r_tol: f64) -> Self {
+    fn r_tol(mut self, r_tol: f64) -> Self {
        self.r_tol = r_tol;
        self
    }
@@ -114,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,
@@ -150,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
@@ -164,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
@@ -183,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
@@ -269,69 +301,94 @@ impl<'a, const D: usize> Integrator<D> for Vern7<D>
 where
    Vern7<D>: Vern7Integrator<'a> + Vern7ExtraStages<'a>,
 {
-    const ORDER: usize = 7;
-    const STAGES: usize = 10;
-    const ADAPTIVE: bool = true;
-    const DENSE: bool = true;
+    fn order(&self) -> usize {
+        7
+    }
+
+    fn stages(&self) -> usize {
+        STAGES
+    }
+
+    fn adaptive(&self) -> bool {
+        true
+    }
+
+    fn dense(&self) -> bool {
+        true
+    }

    // Lazy dense output configuration
-    const MAIN_STAGES: usize = 10;
-    const EXTRA_STAGES: usize = 6;
+    fn main_stages(&self) -> usize {
+        STAGES
+    }

-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
+    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<f64, D>> = vec![SVector::<f64, 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
            let row_start = (i * (i - 1)) / 2;
-            for j in 0..i {
-                y_temp += k[j] * Self::A[row_start + j] * h;
+            for (j, kj) in k.iter().enumerate().take(i) {
+                y_temp += kj * T::lift(Self::A[row_start + j]) * h;
            }
-            k[i] = (ode.f)(ode.t + Self::C[i] * h, y_temp, &ode.params);
+            k[i] = (ode.f)(ode.t + h * Self::C[i], y_temp, &ode.params);
        }

        // Compute 7th order solution using B weights
        let mut next_y = ode.y;
-        for i in 0..Self::STAGES {
-            next_y += k[i] * Self::B[i] * h;
+        for (i, ki) in k.iter().enumerate() {
+            next_y += ki * T::lift(Self::B[i]) * h;
        }

        // Compute error estimate using B_ERROR weights
-        let mut err = SVector::<f64, D>::zeros();
-        for i in 0..Self::STAGES {
-            err += k[i] * Self::B_ERROR[i] * h;
+        let mut err = SVector::<T, D>::zeros();
+        for (i, ki) in k.iter().enumerate() {
+            err += ki * T::lift(Self::B_ERROR[i]) * h;
        }

-        // Compute error norm scaled by tolerance
-        let tol = self.a_tol + ode.y.abs() * self.r_tol;
-        let error_norm = (err.component_div(&tol)).norm();
+        // Compute error norm scaled by tolerance. Step control is a
+        // real-valued decision, so each component drops to the reals as it is
+        // accumulated rather than building real copies of `err` and `y` first.
+        let mut error_norm_sq = 0.0;
+        for i in 0..D {
+            let tol = self.a_tol[i] + ode.y[i].re().abs() * self.r_tol;
+            let scaled = err[i].re() / tol;
+            error_norm_sq += scaled * scaled;
+        }
+        let error_norm = error_norm_sq.sqrt();

        // Store dense output coefficients
        // For now, store all k values for interpolation
        let mut dense_coeffs = vec![ode.y, next_y];
        dense_coeffs.extend_from_slice(&k);

-        (next_y, Some(error_norm), Some(dense_coeffs))
+        StepResult {
+            y: next_y,
+            err: Some(error_norm),
+            dense: Some(dense_coeffs),
+        }
    }

-    fn interpolate(
+    // One interpolation coefficient per line, likewise.
+    #[rustfmt::skip]
+    fn interpolate<T: Scalar>(
        &self,
-        t_start: f64,
-        t_end: f64,
-        dense: &[SVector<f64, D>],
-        t: f64,
-    ) -> SVector<f64, D> {
+        t_start: T,
+        t_end: T,
+        dense: &[SVector<T, D>],
+        t: T,
+    ) -> SVector<T, D> {
        // Vern7 uses 7th order polynomial interpolation
        // Check if extra stages (k11-k16) are available
        // Dense array format: [y0, y1, k1, k2, ..., k10, k11, ..., k16]
@@ -358,8 +415,8 @@ where

        // Helper to evaluate polynomial using Horner's method
        #[inline]
-        fn evalpoly(x: f64, coeffs: &[f64]) -> f64 {
-            let mut result = 0.0;
+        fn evalpoly<T: Scalar>(x: T, coeffs: &[f64]) -> T {
+            let mut result = T::lift(0.0);
            for &c in coeffs.iter().rev() {
                result = result * x + c;
            }
@@ -433,16 +490,16 @@ where
        ]);

        // Compute base interpolation with main stages
-        let mut result = y0 + h * (k1 * b1_theta +
+        let mut result = y0 + (k1 * b1_theta +
                  k4 * b4_theta +
                  k5 * b5_theta +
                  k6 * b6_theta +
                  k7 * b7_theta +
                  k8 * b8_theta +
-                  k9 * b9_theta);
+                  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];
@@ -507,25 +564,25 @@ where
            ]);

            // Add contribution from extra stages
-            result += h * (k11 * b11_theta +
+            result += (k11 * b11_theta +
                          k12 * b12_theta +
                          k13 * b13_theta +
                          k14 * b14_theta +
                          k15 * b15_theta +
-                          k16 * b16_theta);
+                          k16 * b16_theta) * h;
        }

        result
    }

-    fn compute_extra_stages<P>(
+    fn compute_extra_stages<T: Scalar, P>(
        &self,
-        ode: &ODE<D, P>,
-        t_start: f64,
-        y_start: SVector<f64, D>,
-        h: f64,
-        main_stages: &[SVector<f64, D>],
-    ) -> Vec<SVector<f64, D>> {
+        ode: &ODE<T, D, P>,
+        t_start: T,
+        y_start: SVector<T, D>,
+        h: T,
+        main_stages: &[SVector<T, D>],
+    ) -> Vec<SVector<T, D>> {
        // Extract main stages that are used in extra stage computation
        // From Julia: extra stages use k1, k4, k5, k6, k7, k8, k9
        let k1 = &main_stages[0];
@@ -536,90 +593,90 @@ 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;
-        y11 += k1 * Self::A_EXTRA[0] * h;
-        y11 += k4 * Self::A_EXTRA[1] * h;
-        y11 += k5 * Self::A_EXTRA[2] * h;
-        y11 += k6 * Self::A_EXTRA[3] * h;
-        y11 += k7 * Self::A_EXTRA[4] * h;
-        y11 += k8 * Self::A_EXTRA[5] * h;
-        y11 += k9 * Self::A_EXTRA[6] * h;
-        let k11 = (ode.f)(t_start + Self::C_EXTRA[0] * h, y11, &ode.params);
+        y11 += k1 * T::lift(Self::A_EXTRA[0]) * h;
+        y11 += k4 * T::lift(Self::A_EXTRA[1]) * h;
+        y11 += k5 * T::lift(Self::A_EXTRA[2]) * h;
+        y11 += k6 * T::lift(Self::A_EXTRA[3]) * h;
+        y11 += k7 * T::lift(Self::A_EXTRA[4]) * h;
+        y11 += k8 * T::lift(Self::A_EXTRA[5]) * h;
+        y11 += k9 * T::lift(Self::A_EXTRA[6]) * h;
+        let k11 = (ode.f)(t_start + h * Self::C_EXTRA[0], y11, &ode.params);
        extra_stages.push(k11);

        // Stage 12: uses k1, k4-k9, k11 (8 coefficients)
        let mut y12 = y_start;
-        y12 += k1 * Self::A_EXTRA[7] * h;
-        y12 += k4 * Self::A_EXTRA[8] * h;
-        y12 += k5 * Self::A_EXTRA[9] * h;
-        y12 += k6 * Self::A_EXTRA[10] * h;
-        y12 += k7 * Self::A_EXTRA[11] * h;
-        y12 += k8 * Self::A_EXTRA[12] * h;
-        y12 += k9 * Self::A_EXTRA[13] * h;
-        y12 += &extra_stages[0] * Self::A_EXTRA[14] * h; // k11
-        let k12 = (ode.f)(t_start + Self::C_EXTRA[1] * h, y12, &ode.params);
+        y12 += k1 * T::lift(Self::A_EXTRA[7]) * h;
+        y12 += k4 * T::lift(Self::A_EXTRA[8]) * h;
+        y12 += k5 * T::lift(Self::A_EXTRA[9]) * h;
+        y12 += k6 * T::lift(Self::A_EXTRA[10]) * h;
+        y12 += k7 * T::lift(Self::A_EXTRA[11]) * h;
+        y12 += k8 * T::lift(Self::A_EXTRA[12]) * h;
+        y12 += k9 * T::lift(Self::A_EXTRA[13]) * h;
+        y12 += extra_stages[0] * T::lift(Self::A_EXTRA[14]) * h; // k11
+        let k12 = (ode.f)(t_start + h * Self::C_EXTRA[1], y12, &ode.params);
        extra_stages.push(k12);

        // Stage 13: uses k1, k4-k9, k11, k12 (9 coefficients)
        let mut y13 = y_start;
-        y13 += k1 * Self::A_EXTRA[15] * h;
-        y13 += k4 * Self::A_EXTRA[16] * h;
-        y13 += k5 * Self::A_EXTRA[17] * h;
-        y13 += k6 * Self::A_EXTRA[18] * h;
-        y13 += k7 * Self::A_EXTRA[19] * h;
-        y13 += k8 * Self::A_EXTRA[20] * h;
-        y13 += k9 * Self::A_EXTRA[21] * h;
-        y13 += &extra_stages[0] * Self::A_EXTRA[22] * h; // k11
-        y13 += &extra_stages[1] * Self::A_EXTRA[23] * h; // k12
-        let k13 = (ode.f)(t_start + Self::C_EXTRA[2] * h, y13, &ode.params);
+        y13 += k1 * T::lift(Self::A_EXTRA[15]) * h;
+        y13 += k4 * T::lift(Self::A_EXTRA[16]) * h;
+        y13 += k5 * T::lift(Self::A_EXTRA[17]) * h;
+        y13 += k6 * T::lift(Self::A_EXTRA[18]) * h;
+        y13 += k7 * T::lift(Self::A_EXTRA[19]) * h;
+        y13 += k8 * T::lift(Self::A_EXTRA[20]) * h;
+        y13 += k9 * T::lift(Self::A_EXTRA[21]) * h;
+        y13 += extra_stages[0] * T::lift(Self::A_EXTRA[22]) * h; // k11
+        y13 += extra_stages[1] * T::lift(Self::A_EXTRA[23]) * h; // k12
+        let k13 = (ode.f)(t_start + h * Self::C_EXTRA[2], y13, &ode.params);
        extra_stages.push(k13);

        // Stage 14: uses k1, k4-k9, k11, k12, k13 (10 coefficients)
        let mut y14 = y_start;
-        y14 += k1 * Self::A_EXTRA[24] * h;
-        y14 += k4 * Self::A_EXTRA[25] * h;
-        y14 += k5 * Self::A_EXTRA[26] * h;
-        y14 += k6 * Self::A_EXTRA[27] * h;
-        y14 += k7 * Self::A_EXTRA[28] * h;
-        y14 += k8 * Self::A_EXTRA[29] * h;
-        y14 += k9 * Self::A_EXTRA[30] * h;
-        y14 += &extra_stages[0] * Self::A_EXTRA[31] * h; // k11
-        y14 += &extra_stages[1] * Self::A_EXTRA[32] * h; // k12
-        y14 += &extra_stages[2] * Self::A_EXTRA[33] * h; // k13
-        let k14 = (ode.f)(t_start + Self::C_EXTRA[3] * h, y14, &ode.params);
+        y14 += k1 * T::lift(Self::A_EXTRA[24]) * h;
+        y14 += k4 * T::lift(Self::A_EXTRA[25]) * h;
+        y14 += k5 * T::lift(Self::A_EXTRA[26]) * h;
+        y14 += k6 * T::lift(Self::A_EXTRA[27]) * h;
+        y14 += k7 * T::lift(Self::A_EXTRA[28]) * h;
+        y14 += k8 * T::lift(Self::A_EXTRA[29]) * h;
+        y14 += k9 * T::lift(Self::A_EXTRA[30]) * h;
+        y14 += extra_stages[0] * T::lift(Self::A_EXTRA[31]) * h; // k11
+        y14 += extra_stages[1] * T::lift(Self::A_EXTRA[32]) * h; // k12
+        y14 += extra_stages[2] * T::lift(Self::A_EXTRA[33]) * h; // k13
+        let k14 = (ode.f)(t_start + h * Self::C_EXTRA[3], y14, &ode.params);
        extra_stages.push(k14);

        // Stage 15: uses k1, k4-k9, k11, k12, k13 (10 coefficients, reuses k13 not k14)
        let mut y15 = y_start;
-        y15 += k1 * Self::A_EXTRA[34] * h;
-        y15 += k4 * Self::A_EXTRA[35] * h;
-        y15 += k5 * Self::A_EXTRA[36] * h;
-        y15 += k6 * Self::A_EXTRA[37] * h;
-        y15 += k7 * Self::A_EXTRA[38] * h;
-        y15 += k8 * Self::A_EXTRA[39] * h;
-        y15 += k9 * Self::A_EXTRA[40] * h;
-        y15 += &extra_stages[0] * Self::A_EXTRA[41] * h; // k11
-        y15 += &extra_stages[1] * Self::A_EXTRA[42] * h; // k12
-        y15 += &extra_stages[2] * Self::A_EXTRA[43] * h; // k13
-        let k15 = (ode.f)(t_start + Self::C_EXTRA[4] * h, y15, &ode.params);
+        y15 += k1 * T::lift(Self::A_EXTRA[34]) * h;
+        y15 += k4 * T::lift(Self::A_EXTRA[35]) * h;
+        y15 += k5 * T::lift(Self::A_EXTRA[36]) * h;
+        y15 += k6 * T::lift(Self::A_EXTRA[37]) * h;
+        y15 += k7 * T::lift(Self::A_EXTRA[38]) * h;
+        y15 += k8 * T::lift(Self::A_EXTRA[39]) * h;
+        y15 += k9 * T::lift(Self::A_EXTRA[40]) * h;
+        y15 += extra_stages[0] * T::lift(Self::A_EXTRA[41]) * h; // k11
+        y15 += extra_stages[1] * T::lift(Self::A_EXTRA[42]) * h; // k12
+        y15 += extra_stages[2] * T::lift(Self::A_EXTRA[43]) * h; // k13
+        let k15 = (ode.f)(t_start + h * Self::C_EXTRA[4], y15, &ode.params);
        extra_stages.push(k15);

        // Stage 16: uses k1, k4-k9, k11, k12, k13 (10 coefficients, reuses k13 not k14 or k15)
        let mut y16 = y_start;
-        y16 += k1 * Self::A_EXTRA[44] * h;
-        y16 += k4 * Self::A_EXTRA[45] * h;
-        y16 += k5 * Self::A_EXTRA[46] * h;
-        y16 += k6 * Self::A_EXTRA[47] * h;
-        y16 += k7 * Self::A_EXTRA[48] * h;
-        y16 += k8 * Self::A_EXTRA[49] * h;
-        y16 += k9 * Self::A_EXTRA[50] * h;
-        y16 += &extra_stages[0] * Self::A_EXTRA[51] * h; // k11
-        y16 += &extra_stages[1] * Self::A_EXTRA[52] * h; // k12
-        y16 += &extra_stages[2] * Self::A_EXTRA[53] * h; // k13
-        let k16 = (ode.f)(t_start + Self::C_EXTRA[5] * h, y16, &ode.params);
+        y16 += k1 * T::lift(Self::A_EXTRA[44]) * h;
+        y16 += k4 * T::lift(Self::A_EXTRA[45]) * h;
+        y16 += k5 * T::lift(Self::A_EXTRA[46]) * h;
+        y16 += k6 * T::lift(Self::A_EXTRA[47]) * h;
+        y16 += k7 * T::lift(Self::A_EXTRA[48]) * h;
+        y16 += k8 * T::lift(Self::A_EXTRA[49]) * h;
+        y16 += k9 * T::lift(Self::A_EXTRA[50]) * h;
+        y16 += extra_stages[0] * T::lift(Self::A_EXTRA[51]) * h; // k11
+        y16 += extra_stages[1] * T::lift(Self::A_EXTRA[52]) * h; // k12
+        y16 += extra_stages[2] * T::lift(Self::A_EXTRA[53]) * h; // k13
+        let k16 = (ode.f)(t_start + h * Self::C_EXTRA[5], y16, &ode.params);
        extra_stages.push(k16);

        extra_stages
@@ -694,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];
@@ -706,7 +763,8 @@ mod tests {

            while ode.t < t_end {
                let h_step = h.min(t_end - ode.t);
-                let (next_y, _, _) = vern7.step(&ode, h_step);
+                let __step = vern7.step(&ode, h_step);
+                let (next_y, _, _) = (__step.y, __step.err, __step.dense);
                ode.y = next_y;
                ode.t += h_step;
            }
@@ -720,16 +778,8 @@ mod tests {
        let ratio2 = errors[1] / errors[2];

        // Allow some tolerance (expect ratio between 64 and 256)
-        assert!(
-            ratio1 > 64.0 && ratio1 < 256.0,
-            "First ratio: {}",
-            ratio1
-        );
-        assert!(
-            ratio2 > 64.0 && ratio2 < 256.0,
-            "Second ratio: {}",
-            ratio2
-        );
+        assert!(ratio1 > 64.0 && ratio1 < 256.0, "First ratio: {}", ratio1);
+        assert!(ratio2 > 64.0 && ratio2 < 256.0, "Second ratio: {}", ratio2);
    }

    #[test]
@@ -777,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;
@@ -803,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)
@@ -1,20 +1,27 @@
 #![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`.
+#![doc = include_str!("../readme.md")]

 pub mod callback;
 pub mod controller;
 pub mod integrator;
 pub mod ode;
 pub mod problem;
+pub mod scalar;

 pub mod prelude {
-    pub use super::callback::{stop, Callback};
+    pub use super::callback::{stop, Callback, FnCallback};
    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::problem::{Problem, Solution};
+    pub use super::scalar::Scalar;
 }

 #[cfg(test)]
@@ -44,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 = Callback {
-            event: &|t: f64, _y: Vector2<f64>, _p: &Params| 5.0 - t,
+        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
@@ -150,3 +157,241 @@ mod tests {
        );
    }
 }
+
+/// Tests that the derivatives a dual instantiation produces are the right ones.
+#[cfg(test)]
+mod autodiff_tests {
+    use crate::ode::ProblemFunction;
+    use crate::prelude::*;
+    use approx::assert_relative_eq;
+    use nalgebra::{Matrix6, SVector, Vector1, Vector6};
+    use num_dual::{Dual64, DualSVec64};
+    use std::f64::consts::PI;
+
+    const MU: f64 = 3.98600441500000e14;
+    const SMA: f64 = 6.7781363e6;
+
+    /// The two-body right-hand side, shared by the orbit tests below.
+    fn two_body<T: Scalar>(_t: T, state: Vector6<T>, p: &(T,)) -> Vector6<T> {
+        let r = state.fixed_rows::<3>(0);
+        let acc = -(r * p.0) / r.norm().powi(3);
+        Vector6::new(state[3], state[4], state[5], acc[0], acc[1], acc[2])
+    }
+
+    fn orbit_y0<T: Scalar>() -> Vector6<T> {
+        Vector6::new(
+            T::lift(4.263868426884883e6),
+            T::lift(5.146189057155391e6),
+            T::lift(1.1310208421331816e6),
+            T::lift(-5923.454461876975),
+            T::lift(4496.802639690076),
+            T::lift(1870.3893008991558),
+        )
+    }
+
+    fn solve_with<T: Scalar, const D: usize, P, S: Integrator<D> + Copy>(
+        f: ProblemFunction<'_, T, D, P>,
+        integrator: S,
+        t0: T,
+        t_end: T,
+        y0: SVector<T, D>,
+        params: P,
+    ) -> SVector<T, D> {
+        let ode = ODE::new(f, t0, t_end, y0, params);
+        let controller = PIController::new(0.37, 0.04, 10.0, 0.2, 1000.0, 0.9, 0.01);
+        let mut problem = Problem::new(ode, integrator, controller);
+        let solution = problem.solve();
+        *solution.states.last().unwrap()
+    }
+
+    /// The machine-precision case.
+    ///
+    /// A fifth-order Runge-Kutta method integrates `y' = p t^4` exactly, so both
+    /// the state and its derivative with respect to `p` are exact up to
+    /// rounding — there is no truncation error to hide behind. The analytic
+    /// answers are `y(t) = p t^5 / 5` and `dy/dp = t^5 / 5`.
+    #[test]
+    fn test_ad_exact_quadrature() {
+        fn derivative<T: Scalar>(t: T, _y: Vector1<T>, p: &(T,)) -> Vector1<T> {
+            Vector1::new(p.0 * t.powi(4))
+        }
+
+        let t_f = 3.0;
+        // Seed the derivative on the parameter p.
+        let p = Dual64::from_re(2.0).derivative();
+        let y = solve_with(
+            &derivative,
+            DormandPrince45::new().a_tol(1e-12).r_tol(1e-12),
+            Dual64::lift(0.0),
+            Dual64::lift(t_f),
+            Vector1::new(Dual64::lift(0.0)),
+            (p,),
+        );
+
+        assert_relative_eq!(y[0].re, 2.0 * t_f.powi(5) / 5.0, max_relative = 1e-14);
+        assert_relative_eq!(y[0].eps, t_f.powi(5) / 5.0, max_relative = 1e-14);
+    }
+
+    /// `y' = 5y - 3` has the closed-form solution
+    /// `y(t) = (y0 - 0.6) e^{5(t - t0)} + 0.6`, so `dy(t_f)/dy0 = e^{5(t_f - t0)}`
+    /// exactly. Here the comparison is limited by the solver's own truncation
+    /// error rather than by rounding, which is the point: AD returns the exact
+    /// derivative of the numerical solution, and the numerical solution is
+    /// accurate to tolerance.
+    #[test]
+    fn test_ad_matches_analytic_sensitivity() {
+        fn derivative<T: Scalar>(_t: T, y: Vector1<T>, _p: &()) -> Vector1<T> {
+            Vector1::new(y[0] * 5.0 - 3.0)
+        }
+
+        let y0 = Dual64::from_re(1.0).derivative();
+        let y = solve_with(
+            &derivative,
+            DormandPrince45::new().a_tol(1e-14).r_tol(1e-14),
+            Dual64::lift(2.0),
+            Dual64::lift(3.0),
+            Vector1::new(y0),
+            (),
+        );
+
+        assert_relative_eq!(y[0].re, 0.4 * 5.0f64.exp() + 0.6, max_relative = 1e-11);
+        assert_relative_eq!(y[0].eps, 5.0f64.exp(), max_relative = 1e-11);
+    }
+
+    /// A gradient on the orbit problem, checked against a property of the
+    /// two-body flow that is known analytically.
+    ///
+    /// `DualSVec64<6>` carries all six partials at once, so one propagation
+    /// yields the full 6x6 state transition matrix over a period. Two-body
+    /// motion is Hamiltonian, so by Liouville's theorem its flow is
+    /// volume-preserving and `det(Phi) = 1` exactly — for any interval, for any
+    /// mu. A finite-difference Jacobian satisfies this only to the accuracy of
+    /// the difference formula; an exact one satisfies it to the accuracy of the
+    /// propagation.
+    #[test]
+    fn test_ad_orbit_state_transition_matrix() {
+        let period = 2.0 * PI * (SMA.powi(3) / MU).sqrt();
+        let y0: Vector6<f64> = orbit_y0();
+
+        // One propagation in `DualSVec64<6>` yields the whole state transition matrix.
+        let (y_end, stm): (Vector6<f64>, Matrix6<f64>) = num_dual::jacobian(
+            |y: Vector6<DualSVec64<6>>| {
+                type T = DualSVec64<6>;
+                solve_with(
+                    &two_body,
+                    DormandPrince45::new().a_tol(1e-12).r_tol(1e-12),
+                    T::lift(0.0),
+                    T::lift(period),
+                    y,
+                    (T::lift(MU),),
+                )
+            },
+            &y0,
+        );
+
+        // One period returns the state to where it started.
+        for i in 0..6 {
+            assert_relative_eq!(y_end[i], y0[i], max_relative = 1e-10);
+        }
+
+        assert_relative_eq!(stm.determinant(), 1.0, max_relative = 1e-9);
+
+        // And cross-check the whole matrix against central differences, which
+        // agree only to about the square root of machine epsilon.
+        let propagate = |y_start: Vector6<f64>| {
+            solve_with(
+                &two_body,
+                DormandPrince45::new().a_tol(1e-12).r_tol(1e-12),
+                0.0,
+                period,
+                y_start,
+                (MU,),
+            )
+        };
+        for j in 0..6 {
+            let h = 1e-5 * y0[j].abs();
+            let mut plus = y0;
+            let mut minus = y0;
+            plus[j] += h;
+            minus[j] -= h;
+            let fd = (propagate(plus) - propagate(minus)) / (2.0 * h);
+            for i in 0..6 {
+                assert_relative_eq!(stm[(i, j)], fd[i], max_relative = 1e-5, epsilon = 1e-7);
+            }
+        }
+    }
+
+    /// 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.
+    /// That lazy computation calls the right-hand side after the solve has
+    /// finished, so it has to work in dual arithmetic too.
+    #[test]
+    fn test_ad_vern7_lazy_dense_output() {
+        fn derivative<T: Scalar>(_t: T, y: Vector1<T>, _p: &()) -> Vector1<T> {
+            -y
+        }
+
+        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 vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let mut problem = Problem::new(ode, vern7, PIController::default());
+        let solution = problem.solve();
+
+        // Interpolating inside a step is what triggers the extra stages.
+        let t = Dual64::lift(1.3);
+        let y = solution.interpolate(t)[0];
+
+        // y(t) = y0 e^-t, so dy/dy0 = e^-t.
+        assert_relative_eq!(y.re, (-1.3f64).exp(), max_relative = 1e-10);
+        assert_relative_eq!(y.eps, (-1.3f64).exp(), max_relative = 1e-10);
+    }
+}
@@ -1,27 +1,40 @@
 use nalgebra::SVector;

-type ProblemFunction<'a, const D: usize, P> = &'a dyn Fn(f64, SVector<f64, D>, &P) -> SVector<f64, D>;
+use crate::scalar::Scalar;

-/// The basic ODE object that will be passed around. The type (T) and the size (D) will be
+/// The right-hand side of the system: `dy/dt = f(t, y, params)`.
+///
+/// Write it generically over the scalar — `fn f<T: Scalar>(t: T, y: SVector<T, D>, p: &P)` —
+/// and the same function can be instantiated at `f64` to solve the problem and
+/// at a dual number to differentiate it.
+///
+/// Or just write it for f64 if you want to keep things simple
+pub type ProblemFunction<'a, T, const D: usize, P> =
+    &'a dyn Fn(T, SVector<T, D>, &P) -> SVector<T, D>;
+
+/// The basic ODE object that will be passed around. The scalar (T) and the size (D) will be
 /// determined upon creation of the object
 #[derive(Clone, Copy)]
-pub struct ODE<'a, const D: usize, P> {
-    pub f: ProblemFunction<'a, D, P>,
-    pub y: SVector<f64, D>,
-    pub t: f64,
+pub struct ODE<'a, T: Scalar, const D: usize, P> {
+    pub f: ProblemFunction<'a, T, D, P>,
+    pub y: SVector<T, D>,
+    pub t: T,
    pub params: P,
-    pub t0: f64,
-    pub t_end: f64,
-    pub h: f64,
+    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,
 }

-impl<'a, const D: usize, P> ODE<'a, D, P> {
+impl<'a, T: Scalar, const D: usize, P> ODE<'a, T, D, P> {
    pub fn new(
-        f: ProblemFunction<'a, D, P>,
-        t0: f64,
-        t_end: f64,
-        y0: SVector<f64, D>,
+        f: ProblemFunction<'a, T, D, P>,
+        t0: T,
+        t_end: T,
+        y0: SVector<T, D>,
        params: P,
    ) -> Self {
        Self {
@@ -31,10 +44,19 @@ impl<'a, const D: usize, P> ODE<'a, D, P> {
            params,
            t0,
            t_end,
-            h: 0.001,
+            h: T::lift(0.001),
            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)]
@@ -45,7 +67,7 @@ mod tests {
    #[test]
    fn test_ode_creation() {
        type Params = ();
-        fn derivative(_t: f64, y: Vector3<f64>, _p: &Params) -> Vector3<f64> {
+        fn derivative<T: Scalar>(_t: T, y: Vector3<T>, _p: &Params) -> Vector3<T> {
            -y
        }

@@ -64,7 +86,7 @@ mod tests {
        type Params = (f64, bool);
        let params = (34.0, true);

-        fn derivative(t: f64, y: Vector3<f64>, p: &Params) -> Vector3<f64> {
+        fn derivative<T: Scalar>(t: T, y: Vector3<T>, p: &Params) -> Vector3<T> {
            if p.1 {
                -y
            } else {
@@ -0,0 +1,66 @@
+//! The scalar type an ODE is integrated in.
+//!
+//! All of our problem types: state, time and step size, is generic over
+//! [`Scalar`]. Instantiating a problem at `f64` solves it. Instantiating the
+//! same problem at a forward-mode dual number solves it *and* propagates exact
+//! derivatives alongside the state, at machine precision, with no finite
+//! differencing and no extra right-hand-side evaluations beyond the ones the
+//! dual arithmetic performs anyway.
+//!
+//! [`Scalar`] extends the vocabulary [`num_dual`] already defines rather than
+//! 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.
+///
+/// `re`, which discards derivative information, is inherited from
+/// [`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 {
+        Self::from(x)
+    }
+}
+
+impl<T> Scalar for T where
+    T: RealField + Copy + NumOps<f64, T> + DualNum<Primitive = f64> + DualStruct<Real = f64>
+{
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_f64_is_a_scalar() {
+        assert_eq!(f64::lift(3.5), 3.5);
+        assert_eq!(3.5f64.re(), 3.5);
+    }
+
+    #[test]
+    fn test_dual_projects_to_its_real_part() {
+        use num_dual::Dual64;
+
+        // A lifted constant has no derivative; a seeded variable does, and both
+        // project back to the same real value.
+        let c = <Dual64 as Scalar>::lift(3.5);
+        assert_eq!(c.re, 3.5);
+        assert_eq!(c.eps, 0.0);
+
+        let x = Dual64::from_re(3.5).derivative();
+        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);
+    }
+}
Author	SHA1	Message	Date
connor	981b1c8da9	Merge branch 'release/0.8.0'	2026-10-02 19:44:57 -04:00
connor	95f19f1a48	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().	2026-10-02 19:44:52 -04:00
connor	e25d654ef1	Simplified scalar type to just use num_dual	2026-10-02 19:44:00 -04:00
connor	6f215788e5	Merge branch 'release/0.7.0'	2026-10-01 14:14:12 -04:00
connor	17105964a2	Release 0.7.0 A Problem can own a callback as well as borrow one. with_owned_callback takes a Box<dyn Callback + 'o>, which is what a type that stores the Problem needs: it cannot lend a callback from one of its own fields without becoming self-referential. Problem gains a second lifetime for this, Problem<'a, 'o, T, D, S, P, C>, so code that names the type explicitly needs Problem<'a, '_, ..>. Both are inferred and no call site changes. The lifetimes have to be separate because dropping a Box<dyn Callback> runs a destructor the compiler cannot see into, and sharing one lifetime would keep a borrowed callback borrowed until the problem went out of scope, which is exactly what the borrowed form exists to avoid. Event location, root finding, jump recording and the controller reset are unchanged.	2026-10-01 14:13:33 -04:00
connor	0e4a0b1798	Added an ability for an owned callback	2026-10-01 14:12:51 -04:00
connor	78ac6a60a0	Merge branch 'release/0.6.1'	2026-10-01 13:09:00 -04:00
connor	68ddc431a5	Release 0.6.1 Documentation and housekeeping only; the library behaves exactly as 0.6.0 does. Three rustdoc warnings in the Jacobian comment are gone, so the formulas there render as code rather than as attempted links, the tree is rustfmt-clean with the Butcher tableaux and interpolation coefficients marked to keep their hand alignment, and the git hooks stay out of the package.	2026-10-01 13:08:11 -04:00
connor	144e74c517	Merge branch 'fix/package-excludes'	2026-10-01 13:07:58 -04:00
connor	a1fbcf33cf	Keep the git hooks out of the published package exclude covered roadmap/ but not hooks/, so the pre-commit and pre-push scripts would have gone out in the next crate tarball. They are for working in the repo and mean nothing to a consumer of the library.	2026-10-01 13:07:17 -04:00
connor	cd3c40940e	Git hooks, format changes, etc.	2026-10-01 13:05:59 -04:00
connor	328abfde1a	Merge branch 'release/0.6.0'	2026-10-01 12:47:36 -04:00
connor	a73e3dc9ec	Release 0.6.0 Callback effects can change the state. The effect runs after the step has been shortened to the crossing and committed, so what it leaves in ode.y is what the integration continues from; it was overwritten immediately before. Callback::effect takes &mut self, Problem::with_callback takes &mut dyn Callback, and Controller requires a reset method, all of which are breaking. A jump records the state on both sides of itself and resets the step size controller, interpolating at one reports the state the effect left, stop ends the solution at the crossing without disturbing t_end, and the earliest of two events in one step is the one that fires.	2026-10-01 12:47:33 -04:00
connor	c5768ad9c0	Fixed event handling to properly keep the state changes, plus a few other fixes	2026-10-01 12:41:39 -04:00
connor	55884190fa	Merge branch 'release/0.5.0'	2026-09-30 14:08:30 -04:00
connor	0b6996b0fc	Release 0.5.0 Runtime integrator selection. The structural facts on Integrator move from associated constants to &self methods, which is a breaking change for anyone implementing the trait or reading ORDER, STAGES, ADAPTIVE or DENSE. In exchange Method dispatches to any of the four integrators as a value, and MethodName names one in a configuration file.	2026-09-30 14:08:28 -04:00
connor	5b40932317	Merge branch 'feature/runtime-integrator-selection'	2026-09-30 14:08:16 -04:00
connor	8894de37d6	Added the ability to select a method at runtime if so desired	2026-09-30 14:06:49 -04:00
connor	2ce5028529	Merge branch 'release/0.4.1'	2026-09-30 13:49:43 -04:00
connor	3e6f79b9bd	Release 0.4.1 The Tolerances trait: the a_tol, a_tol_full and r_tol setters move off the inherent impls onto a trait, with absolute_tolerance and relative_tolerance getters, so code generic over the integrator can set and read them.	2026-09-30 13:49:34 -04:00
connor	bbc94592e5	Merge branch 'feature/tolerances-trait'	2026-09-30 13:49:14 -04:00
connor	da479a378e	Made the tolerance trait separable from the integrator itself	2026-09-30 13:48:20 -04:00
connor	812b5adaf9	Merge branch 'release/0.4.0'	2026-09-29 18:29:35 -04:00
connor	8155a561ac	Release 0.4.0 Generic scalar type and forward-mode autodiff, the step size controller fix, and the clippy and benchmark-report work that came with them.	2026-09-29 18:29:33 -04:00
connor	788f93cccb	Merge pull request 'Genericization and Autodiff' (#4 ) from feature/generic-scalar into main Reviewed-on: #4	2026-09-29 18:24:02 -04:00
connor	2053aaa68e	Updated reports after bug fix	2026-09-29 18:21:40 -04:00
connor	c20a090036	Fixed a small bug on the step sizing	2026-09-29 18:02:58 -04:00
connor	db14a02df4	Updated docs on benchmarks	2026-09-29 18:00:09 -04:00
connor	f8a47d0126	Incorporated the other features: genericization and autodiff mainly	2026-09-29 17:46:33 -04:00
connor	cbbe455c07	Merge branch 'fix/release-metadata'	2026-09-29 16:08:05 -04:00
connor	1f2070d110	Fix release metadata and repair the benchmarks - repository pointed at the old host and the old repo name; this is the link crates.io and docs.rs display - the readme example imported `differential_equations`, so it could not compile under the crate's published name - `cargo bench` did not build: `Solution` borrows the ODE, so `.solve()` cannot be called on a temporary `Problem`. Bind it first.	2026-09-29 16:08:03 -04:00