Made the tolerance trait separable from the integrator itself

Merge branch 'release/0.4.0'
Release 0.4.0
2026-09-30 13:48:20 -04:00 · 2026-09-29 18:29:35 -04:00 · 2026-09-29 18:29:33 -04:00 · 2026-09-29 18:24:02 -04:00 · 2026-09-29 18:21:40 -04:00 · 2026-09-29 18:02:58 -04:00
24 changed files with 1870 additions and 664 deletions
@@ -1,24 +1,28 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.3.0"
+version = "0.4.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
+[features]
+default = ["autodiff"]
+autodiff = ["dep:num-dual"]

 [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 = { version = "0.15", optional = true }
 roots = "0.0.8"

 [dev-dependencies]
 approx = "0.5"
+num-dual = "0.15"
 criterion = "0.7.0"

 [[bench]]
@@ -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,10 @@ 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 +203,10 @@ 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 +233,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 +242,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());
            });
        });
    }
@@ -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());
        });
    });
 }
@@ -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,23 @@ 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 +211,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 +220,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());
            });
        });
    }
@@ -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,8 +15,8 @@ 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.

@@ -27,6 +27,8 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy
 - **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

 ### Future Improvements

@@ -47,7 +49,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 +70,77 @@ 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 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(&value_too_high);
 let solution = problem.solve();

 // Can interpolate solutions to whatever you want
 let _interpolated_answer = solution.interpolate(4.4);
 ```
+
+## Automatic differentiation
+
+State, time and step size are generic over the `Scalar` trait. `f64` implements it by default,
+and so do the forward-mode dual numbers of [`num-dual`](https://crates.io/crates/num-dual) (the
+`autodiff` feature, on by default). 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.
@@ -58,6 +58,23 @@ 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
+
+- [ ] **[Runtime Integrator Selection](features/41-runtime-integrator-selection.md)**
+  - 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
+
 ### Callbacks

 - [ ] **[Discrete Callbacks](features/05-discrete-callbacks.md)**
@@ -298,6 +315,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 +353,18 @@ 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: 41
+- Tier 1: 10 features (5/10 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: 12.2% (5/41 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)

-Last updated: 2025-10-24
+Last updated: 2026-09-30
@@ -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,174 @@
+# Feature: Runtime Integrator Selection
+
+**Status**: proposed (2026-09-30)
+
+## 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
@@ -1,21 +1,54 @@
 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 callback function you can utilize to change your model mid-run without stopping.
 ///
-/// The integration solver will check this function for zero crossings
+/// The solver evaluates [`event`](Callback::event) at each end of every accepted
+/// step and looks for a sign change. When it finds one it locates the crossing
+/// and applies [`effect`](Callback::effect) to the ODE there.
+/// ```
+/// # use ordinary_diffeq::prelude::*;
+/// # use nalgebra::Vector2;
+/// # type Params = ();
+/// let too_high = FnCallback {
+///     event: &|_t: f64, y: &Vector2<f64>, _p: &Params| 10.0 - y[0],
+///     effect: &stop,
+/// };
+/// ```
+pub trait Callback<T: Scalar, const D: usize, P> {
+    /// The function whose zero crossings mark the event.
+    fn event(&self, t: T, y: &SVector<T, D>, params: &P) -> T;
+
+    /// What to do to the ODE once the crossing has been located.
+    fn effect(&self, ode: &mut ODE<'_, T, D, P>);
+}
+
+/// A [`Callback`] built from a pair of closures.
 #[derive(Clone, Copy)]
-pub struct Callback<'a, const D: usize, P> {
+pub struct FnCallback<'a, T: Scalar, const D: usize, P> {
    /// The function to check for zero crossings
-    pub event: &'a dyn Fn(f64, SVector<f64, D>, &P) -> f64,
+    pub event: &'a dyn Fn(T, &SVector<T, D>, &P) -> T,

    /// The function to change the ODE
-    pub effect: &'a dyn Fn(&mut ODE<D, P>),
+    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(&self, ode: &mut ODE<'_, T, D, P>) {
+        (self.effect)(ode)
+    }
 }

 /// A convenience function for stopping the integration
-pub fn stop<const D: usize, P>(ode: &mut ODE<D, P>) {
+pub fn stop<T: Scalar, const D: usize, P>(ode: &mut ODE<'_, T, D, P>) {
    ode.t_end = ode.t;
 }

@@ -26,8 +59,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,18 @@ 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>);
 }

 #[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 +45,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 +61,30 @@ 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;
+    }
+
 }

-impl PIController {
+impl<T: Scalar> PIController<T> {
    pub fn new(
        alpha: f64,
        beta: f64,
@@ -83,17 +102,18 @@ 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)
    }
 }

+
 /// PID (Proportional-Integral-Derivative) step size controller
 ///
 /// The PID controller is an advanced adaptive time-stepping controller that provides
@@ -111,9 +131,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,11 +158,11 @@ 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)
@@ -156,35 +176,35 @@ pub struct PIDController {
    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);
+            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)
            )
        } 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
            };
@@ -200,43 +220,53 @@ 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.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;
+    }
+
 }

-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,9 +274,9 @@ 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(),
@@ -254,16 +284,16 @@ impl PIDController {
            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),
+            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,9 +301,9 @@ 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
+            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
@@ -292,22 +322,22 @@ impl PIDController {
    }
 }

-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)
+            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
@@ -317,13 +347,14 @@ impl Default for PIDController {
    }
 }

+
 #[cfg(test)]
 mod tests {
    use super::*;

    #[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 +368,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 +384,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 +424,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 +435,8 @@ 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 +469,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 +495,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 +518,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 +532,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 +547,8 @@ 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();
@@ -550,7 +582,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,8 @@
 use nalgebra::SVector;

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

 /// Bogacki-Shampine 3/2 integrator trait for tableau coefficients
 pub trait BS3Integrator<'a> {
@@ -57,6 +58,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 +80,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
    }
@@ -140,38 +160,52 @@ where
    const ADAPTIVE: bool = true;
    const DENSE: bool = true;

-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
+    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(); Self::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
@@ -183,16 +217,20 @@ where
            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 +253,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)
    }
 }

@@ -233,6 +271,7 @@ mod tests {
    use nalgebra::Vector1;

    #[test]
+    #[allow(clippy::assertions_on_constants)]
    fn test_bs3_creation() {
        let _bs3: BS3<1> = BS3::new();
        assert_eq!(BS3::<1>::ORDER, 3);
@@ -254,7 +293,8 @@ mod tests {
        let bs3 = BS3::new();
        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();
@@ -283,7 +323,8 @@ mod tests {
        let bs3 = BS3::new();
        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 +353,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 +387,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 +432,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 +441,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,8 @@
 use nalgebra::SVector;

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

 /// Integrator Trait
 pub trait DormandPrinceIntegrator<'a> {
@@ -11,9 +12,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 +68,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()
    }
 }

@@ -96,51 +179,58 @@ where
    const ADAPTIVE: bool = true;
    const DENSE: bool = true;

-    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 step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
+        let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
        let mut 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[Self::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 {
            // 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 + Self::STAGES]) * h;
+            rcont5 += k[i] * h * T::lift(Self::D[i]);
        }
        let rcont2 = next_y - ode.y;
-        let rcont3 = h * k[0] - rcont2;
+        let rcont3 = k[0] * h - rcont2;
        let rcont4 = rcont2 - k[Self::STAGES - 1] * h - rcont3;
-        let tol = self.a_tol + ode.y * self.r_tol;
        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
    }
 }
@@ -1,12 +1,34 @@
 use nalgebra::SVector;

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

 pub mod bs3;
 pub mod dormand_prince;
 pub mod rosenbrock;
 pub mod vern7;

+/// 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>,
+
+    /// 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
 pub trait Integrator<const D: usize> {
    const ORDER: usize;
@@ -23,21 +45,18 @@ pub trait Integrator<const D: usize> {
    /// Total stages when full dense output is computed
    const TOTAL_DENSE_STAGES: 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 +75,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 +138,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 +151,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)]
@@ -31,24 +32,26 @@ impl Default for JacobianUpdateStrategy {
 ///
 /// 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>,
+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 +134,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 +143,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 +165,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 +175,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 +197,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()
@@ -218,81 +234,60 @@ impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
    const ADAPTIVE: bool = true;
    const DENSE: bool = true;

-    fn step<P>(
-        &self,
-        ode: &ODE<D, P>,
-        h: f64,
-    ) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
+    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 +295,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 +322,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 +382,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 +416,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;
            }
@@ -465,7 +466,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;
            }
@@ -497,7 +499,6 @@ mod tests {

    #[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;
@@ -518,7 +519,6 @@ mod tests {

        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());
@@ -1,7 +1,8 @@
 use nalgebra::SVector;

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

 /// Verner 7 integrator trait for tableau coefficients
 pub trait Vern7Integrator<'a> {
@@ -69,6 +70,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 +92,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
    }
@@ -278,13 +298,9 @@ where
    const MAIN_STAGES: usize = 10;
    const EXTRA_STAGES: usize = 6;

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

        // Stage 1: k[0] = f(t, y)
        k[0] = (ode.f)(ode.t, ode.y, &ode.params);
@@ -295,43 +311,54 @@ where
            // 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(
+    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 +385,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,13 +460,13 @@ 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 {
@@ -507,25 +534,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];
@@ -540,86 +567,86 @@ where

        // 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
@@ -706,7 +733,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;
            }
@@ -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`. Requires `autodiff`
+// because one of those examples differentiates.
+#![cfg_attr(feature = "autodiff", 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::rosenbrock::Rosenbrock23;
    pub use super::integrator::vern7::Vern7;
    pub use super::ode::ODE;
+    pub use super::integrator::{Integrator, StepResult, Tolerances};
    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 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(&value_too_high);
        let solution = problem.solve();

        // Can interpolate solutions to whatever you want
@@ -150,3 +157,193 @@ mod tests {
        );
    }
 }
+
+/// Tests that the derivatives a dual instantiation produces are the right ones.
+#[cfg(all(test, feature = "autodiff"))]
+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);
+            }
+        }
+    }
+
+    /// 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,37 @@
 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,
    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,7 +41,7 @@ impl<'a, const D: usize, P> ODE<'a, D, P> {
            params,
            t0,
            t_end,
-            h: 0.001,
+            h: T::lift(0.001),
            finished: false,
        }
    }
@@ -45,7 +55,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 +74,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 {
@@ -1,28 +1,33 @@
+//! Assembling an ODE, an integrator and a controller
+
 use nalgebra::SVector;
 use roots::{find_root_brent, SimpleConvergency};
 use std::cell::RefCell;

 use super::callback::Callback;
-use super::controller::{Controller, PIController, TryStep};
+use super::controller::{Controller, TryStep};
 use super::integrator::Integrator;
 use super::ode::ODE;
+use super::scalar::Scalar;

 #[derive(Clone)]
-pub struct Problem<'a, const D: usize, S, P>
+pub struct Problem<'a, T: Scalar, const D: usize, S, P, C>
 where
    S: Integrator<D>,
+    C: Controller<T>,
 {
-    ode: ODE<'a, D, P>,
+    ode: ODE<'a, T, D, P>,
    integrator: S,
-    controller: PIController,
-    callbacks: Vec<Callback<'a, D, P>>,
+    controller: C,
+    callbacks: Vec<&'a dyn Callback<T, D, P>>,
 }

-impl<'a, const D: usize, S, P> Problem<'a, D, S, P>
+impl<'a, T: Scalar, const D: usize, S, P, C> Problem<'a, T, D, S, P, C>
 where
    S: Integrator<D> + Copy,
+    C: Controller<T>,
 {
-    pub fn new(ode: ODE<'a, D, P>, integrator: S, controller: PIController) -> Self {
+    pub fn new(ode: ODE<'a, T, D, P>, integrator: S, controller: C) -> Self {
        Problem {
            ode,
            integrator,
@@ -30,69 +35,79 @@ where
            callbacks: Vec::new(),
        }
    }
-    pub fn solve(&mut self) -> Solution<'_, S, D, P> {
+
+    pub fn solve(&mut self) -> Solution<'_, T, S, D, P> {
        let mut convergency = SimpleConvergency {
-            eps: 1e-12,
+            eps: 1e-12f64,
            max_iter: 1000,
        };
-        let mut times: Vec<f64> = vec![self.ode.t];
-        let mut states: Vec<SVector<f64, D>> = vec![self.ode.y];
-        let mut dense_coefficients: Vec<RefCell<Vec<SVector<f64, D>>>> = Vec::new();
-        while self.ode.t < self.ode.t_end {
-            if self.ode.t + self.controller.next_step_guess.extract() > self.ode.t_end {
+        let mut times: Vec<T> = vec![self.ode.t];
+        let mut states: Vec<SVector<T, D>> = vec![self.ode.y];
+        let mut dense_coefficients: Vec<RefCell<Vec<SVector<T, D>>>> = Vec::new();
+        while self.ode.t.re() < self.ode.t_end.re() {
+            if (self.ode.t + self.controller.next_step_guess().extract()).re()
+                > self.ode.t_end.re()
+            {
                // If the next step would go past the end, then just set it to the end
-                self.controller.next_step_guess = TryStep::NotYetAccepted(
-                    self.ode.t_end - self.ode.t,
-                );
+                self.controller
+                    .set_next_step_guess(TryStep::NotYetAccepted(self.ode.t_end - self.ode.t));
            }
            let (mut new_y, mut curr_step, mut dense_option) = if S::ADAPTIVE {
                // First, we try stepping with the "next step guess" to get the error
-                let (mut trial_y, mut err_option, mut dense_option) =
-                    self.integrator.step(&self.ode, self.controller.next_step_guess.extract());
-                let mut err = err_option.unwrap();
+                let mut result = self
+                    .integrator
+                    .step(&self.ode, self.controller.next_step_guess().extract());
+                let mut err = result.err.unwrap();
                // Then we determine whether we need to reduce the step size or not
                // If successful, we get the next step guess
-                let initial_guess = self.controller.next_step_guess.extract();
-                let mut next_step_guess = <PIController as Controller<D>>::determine_step(
-                    &mut self.controller,
-                    initial_guess,
-                    err,
-                );
+                let initial_guess = self.controller.next_step_guess().extract();
+                let mut next_step_guess = self.controller.determine_step(initial_guess, err);
                while !next_step_guess.is_accepted() {
                    // If that step isn't acceptable, then change the step until it is
-                    (trial_y, err_option, dense_option) =
-                        self.integrator.step(&self.ode, next_step_guess.extract());
-                    next_step_guess = <PIController as Controller<D>>::determine_step(
-                        &mut self.controller,
-                        next_step_guess.extract(),
-                        err,
-                    );
-                    err = err_option.unwrap();
+                    result = self.integrator.step(&self.ode, next_step_guess.extract());
+                    err = result.err.unwrap();
+                    next_step_guess = self
+                        .controller
+                        .determine_step(next_step_guess.extract(), err);
                }
                // So at this point we can safely assume we have an accepted step
-                self.controller.next_step_guess = next_step_guess.reset().unwrap();
-                (trial_y, next_step_guess.extract(), dense_option)
+                self.controller
+                    .set_next_step_guess(next_step_guess.reset().unwrap());
+                (result.y, next_step_guess.extract(), result.dense)
            } else {
                // If fixed time step just step forward one step
-                let (trial_y, _, dense_option) = self.integrator.step(&self.ode, self.controller.next_step_guess.extract());
-                (trial_y, self.controller.next_step_guess.extract(), dense_option)
+                let result = self
+                    .integrator
+                    .step(&self.ode, self.controller.next_step_guess().extract());
+                (
+                    result.y,
+                    self.controller.next_step_guess().extract(),
+                    result.dense,
+                )
            };
            if !self.callbacks.is_empty() {
                // Check for events occurring
                for callback in &self.callbacks {
-                    if (callback.event)(self.ode.t, self.ode.y, &self.ode.params)
-                        * (callback.event)(self.ode.t + curr_step, new_y, &self.ode.params)
-                        < 0.0
-                    {
+                    // The event condition is a real-valued question, so it is
+                    // asked in the reals — see this module's header.
+                    let before = callback.event(self.ode.t, &self.ode.y, &self.ode.params).re();
+                    let after = callback
+                        .event(self.ode.t + curr_step, &new_y, &self.ode.params)
+                        .re();
+                    if before * after < 0.0 {
                        // If the event crossed zero, then find the root
-                        let f = |test_t| {
-                            let test_y = self.integrator.step(&self.ode, test_t).0;
-                            (callback.event)(self.ode.t + test_t, test_y, &self.ode.params)
+                        let f = |test_h: f64| {
+                            let h = T::lift(test_h);
+                            let test_y = self.integrator.step(&self.ode, h).y;
+                            callback.event(self.ode.t + h, &test_y, &self.ode.params).re()
                        };
-                        let root = find_root_brent(0.0, curr_step, &f, &mut convergency).unwrap();
-                        curr_step = root;
-                        (new_y, _, dense_option) = self.integrator.step(&self.ode, curr_step);
-                        (callback.effect)(&mut self.ode);
+                        let root =
+                            find_root_brent(0.0f64, curr_step.re(), &f, &mut convergency).unwrap();
+                        curr_step = T::lift(root);
+                        let result = self.integrator.step(&self.ode, curr_step);
+                        new_y = result.y;
+                        dense_option = result.dense;
+                        callback.effect(&mut self.ode);
                    }
                }
            }
@@ -112,52 +127,48 @@ where
        }
    }

-    pub fn with_callback(mut self, callback: Callback<'a, D, P>) -> Self {
+    pub fn with_callback(mut self, callback: &'a dyn Callback<T, D, P>) -> Self {
        self.callbacks.push(callback);
-        Self {
-            ode: self.ode,
-            integrator: self.integrator,
-            controller: self.controller,
-            callbacks: self.callbacks,
-        }
+        self
    }
 }

-pub struct Solution<'a, S, const D: usize, P>
+pub struct Solution<'a, T: Scalar, S, const D: usize, P>
 where
    S: Integrator<D>,
 {
-    pub ode: &'a ODE<'a, D, P>,
+    pub ode: &'a ODE<'a, T, D, P>,
    pub integrator: S,
-    pub times: Vec<f64>,
-    pub states: Vec<SVector<f64, D>>,
-    pub dense: Vec<RefCell<Vec<SVector<f64, D>>>>,
+    pub times: Vec<T>,
+    pub states: Vec<SVector<T, D>>,
+    pub dense: Vec<RefCell<Vec<SVector<T, D>>>>,
 }

-impl<'a, S, const D: usize, P> Solution<'a, S, D, P>
+impl<T: Scalar, S, const D: usize, P> Solution<'_, T, S, D, P>
 where
    S: Integrator<D>,
 {
-    pub fn interpolate(&self, t: f64) -> SVector<f64, D> {
+    pub fn interpolate(&self, t: T) -> SVector<T, D> {
        // First check that the t is within bounds
-        let last = self.times.last().unwrap();
-        let first = self.times.first().unwrap();
+        let last = self.times.last().unwrap().re();
+        let first = self.times.first().unwrap().re();

        // TODO: Improve these errors
-        let mut times = self.times.clone();
-        if *first > *last {
+        let mut times: Vec<f64> = self.times.iter().map(|x| x.re()).collect();
+        if first > last {
            times.reverse();
        }
-        if t < *first || t > *last {
+        if t.re() < first || t.re() > last {
            panic!();
        }

        // Then find the two t values closest to the desired t
-        match times.binary_search_by(|x| x.total_cmp(&t)) {
+        match times.binary_search_by(|x| x.total_cmp(&t.re())) {
            Ok(index) => self.states[index],
            Err(end_index) => {
-                let t_start = times[end_index - 1];
-                let t_end = times[end_index];
+                // Then send that to the integrator
+                let t_start = self.times[end_index - 1];
+                let t_end = self.times[end_index];
                let y_start = self.states[end_index - 1];
                let h = t_end - t_start;

@@ -205,9 +216,10 @@ where
 #[cfg(test)]
 mod tests {
    use super::*;
-    use crate::callback::stop;
+    use crate::callback::{stop, FnCallback};
    use crate::controller::PIController;
    use crate::integrator::dormand_prince::DormandPrince45;
+    use crate::integrator::Tolerances;
    use approx::assert_relative_eq;
    use nalgebra::Vector3;

@@ -247,12 +259,12 @@ mod tests {
        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-5);
        let controller = PIController::default();

-        let value_too_high = Callback {
-            event: &|_: f64, y: SVector<f64, 3>, _: &Params| 10.0 - y[0],
+        let value_too_high = FnCallback {
+            event: &|_: f64, y: &SVector<f64, 3>, _: &Params| 10.0 - y[0],
            effect: &stop,
        };

-        let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high);
+        let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high);
        let solution = problem.solve();

        assert_relative_eq!(
@@ -0,0 +1,119 @@
+//! 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`] is implemented for `f64` always, and for the dual numbers of the
+//! [`num_dual`] crate under the `autodiff` feature (on by default).
+
+use nalgebra::RealField;
+use num_traits::NumOps;
+
+/// A scalar that an ODE's state, time and step size may take.
+pub trait Scalar: RealField + Copy + NumOps<f64, Self> {
+    /// This value with all derivative information discarded.
+    fn re(&self) -> f64;
+
+    /// A real constant as a value of this type, with zero derivative.
+    fn lift(x: f64) -> Self;
+}
+
+impl Scalar for f64 {
+    #[inline(always)]
+    fn re(&self) -> f64 {
+        *self
+    }
+
+    #[inline(always)]
+    fn lift(x: f64) -> Self {
+        x
+    }
+}
+
+#[cfg(feature = "autodiff")]
+mod dual {
+    use super::Scalar;
+    use nalgebra::{Const, RealField};
+    use num_dual::{DualNum, DualStruct};
+
+    /// Implements `Scalar` for one of `num-dual`'s dual number types.
+    macro_rules! impl_scalar {
+        ($ty:ident) => {
+            impl<T> Scalar for num_dual::$ty<T>
+            where
+                Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
+            {
+                #[inline(always)]
+                fn re(&self) -> f64 {
+                    DualStruct::re(self)
+                }
+
+                #[inline(always)]
+                fn lift(x: f64) -> Self {
+                    Self::from(x)
+                }
+            }
+        };
+        ($ty:ident, $($n:ident),+) => {
+            impl<T: nalgebra::Scalar, $(const $n: usize),+> Scalar
+                for num_dual::$ty<T, $(Const<$n>),+>
+            where
+                Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
+            {
+                #[inline(always)]
+                fn re(&self) -> f64 {
+                    DualStruct::re(self)
+                }
+
+                #[inline(always)]
+                fn lift(x: f64) -> Self {
+                    Self::from(x)
+                }
+            }
+        };
+    }
+
+    impl_scalar!(Dual);
+    impl_scalar!(Dual2);
+    impl_scalar!(Dual3);
+    impl_scalar!(HyperDual);
+    impl_scalar!(HyperHyperDual);
+    impl_scalar!(DualVec, N);
+    impl_scalar!(Dual2Vec, N);
+    impl_scalar!(HyperDualVec, M, N);
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_f64_is_a_scalar() {
+        assert_eq!(f64::lift(3.5), 3.5);
+        assert_eq!(Scalar::re(&3.5), 3.5);
+    }
+
+    #[cfg(feature = "autodiff")]
+    #[test]
+    fn test_dual_projects_to_its_real_part() {
+        use num_dual::{Dual64, DualNum};
+
+        // 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!(Scalar::re(&x), 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	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