Merge branch 'release/0.8.0'

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

num-dual is a plain dependency rather than optional.

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

Cargo.toml

@@ -1,28 +1,23 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.6.0"
+version = "0.8.0"
 edition = "2021"
 authors = ["Connor Johnstone"]
 description = "A library for solving differential equations based on the DifferentialEquations.jl julia library."
 readme = "readme.md"
 repository = "https://git.rcjohnstone.com/connor/ordinary-diffeq"
 license = "MIT"
-exclude = ["roadmap/"]
-
-[features]
-default = ["autodiff"]
-autodiff = ["dep:num-dual"]
+exclude = ["roadmap/", "hooks/"]
 
 [dependencies]
 serde = { version = "1.0", features = ["derive"] }
 nalgebra = { version = "0.35", features = ["serde-serialize"] }
 num-traits = "0.2.19"
-num-dual = { version = "0.15", optional = true }
+num-dual = "0.15"
 roots = "0.0.8"
 
 [dev-dependencies]
 approx = "0.5"
-num-dual = "0.15"
 criterion = "0.7.0"
 serde_json = "1.0"
 

benches/bs3_vs_dp5.rs

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

benches/orbit.rs

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

benches/simple_1d.rs

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

benches/vern7_comparison.rs

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

hooks/pre-commit

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

hooks/pre-push

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

readme.md

@@ -134,6 +134,58 @@ crossing as its last point and the problem's own `t_end` left alone.
 If two events fall inside one step the earlier one fires and the step ends there; the later one is
 found again on the step after it.
 
+### Borrowed or owned
+
+`with_callback` borrows, which is what lets you read a stateful callback's own fields once the solve
+is done:
+
+```rust
+# use nalgebra::Vector1;
+# use ordinary_diffeq::prelude::*;
+struct Count {
+    fired: usize,
+}
+
+impl Callback<f64, 1, ()> for Count {
+    fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
+        1.0 - t
+    }
+
+    fn effect(&mut self, _ode: &mut ODE<'_, f64, 1, ()>) {
+        self.fired += 1;
+    }
+}
+
+fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+    Vector1::new(1.0)
+}
+
+let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(0.0), ());
+let mut counter = Count { fired: 0 };
+let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+    .with_callback(&mut counter);
+problem.solve();
+
+assert_eq!(counter.fired, 1);
+```
+
+`with_owned_callback` takes a `Box` instead, for a caller that has nothing to lend. A type that wraps
+this one has to keep the `Problem` in a field, since `solve` hands back a solution borrowing from it,
+and it then cannot also hold a callback in a field and lend it to the problem — that is a
+self-referential struct, and moving the wrapper would dangle the reference. Moving the callback into
+the problem sidesteps it:
+
+```rust,ignore
+let problem = Problem::new(ode, integrator, controller)
+    .with_owned_callback(Box::new(my_callback));
+self.problem = Some(problem);
+self.problem.as_mut().unwrap().solve()
+```
+
+An owned callback may hold borrows of its own, so a boxed `FnCallback` works. Nothing can read its
+state afterwards; a callback with something to report either stays borrowed, or shares an
+`Rc<Cell<_>>` with its owner, which can be read while the solution is still alive.
+
 ## Choosing the integrator at runtime
 
 The integrator is normally a type parameter, fixed when the code is compiled. `Method` is an enum
@@ -166,10 +218,13 @@ holds several methods for comparison over the same problem.
 
 ## Automatic differentiation
 
-State, time and step size are generic over the `Scalar` trait. `f64` implements it by default,
-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.
+State, time and step size are generic over the `Scalar` trait, which is a blanket extension of
+`DualNum` and `DualStruct` from [`num-dual`](https://crates.io/crates/num-dual): those two supply
+the arithmetic, the transcendentals and `re`, which discards derivative information, and `Scalar`
+adds only `lift`. Both are implemented for `f64` as well as for every one of that crate's
+forward-mode dual numbers, so `Scalar` covers all of them, and a crate depending on this one
+alongside another built on `num-dual` can satisfy both with one bound. Solving at a dual scalar
+propagates exact derivatives alongside the state.
 
 The only change needed is to the signature of the right-hand side. Here is the pendulum from above,
 before and after:

roadmap/GETTING_STARTED.md

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

roadmap/README.md

@@ -86,6 +86,8 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - `effect(&mut self)`, termination at the crossing, both sides of a jump
     recorded, the controller reset across one, and a defined order for two events
     in one step
+  - extended in 0.7.0 so a `Problem` can own a callback as well as borrow one,
+    which is what a type that stores the problem needs
   - Raised by `astrokit`, which needs a Modified Rodrigues Parameter shadow set
     swap mid-propagation
   - **Dependencies**: None

roadmap/features/42-callback-effects.md

@@ -1,6 +1,6 @@
 # Feature: Callback Effects That Can Change the State
 
-**Status**: completed (2026-10-01)
+**Status**: completed (2026-10-01), extended with owned callbacks (2026-10-01)
 
 ## Overview
 
@@ -166,6 +166,108 @@ Brent bracketing a root at `h = 0`; and refusing a located crossing below some
 `h`, which is the same arbitrary fraction one step further down. If something
 cleaner turns up, this is the place to change it.
 
+## Follow-up: a problem can own a callback (0.7.0)
+
+`Problem` held `Vec<&'a mut dyn Callback>`, so something else had to own each
+callback and outlive the problem. That is fine from a function body and impossible
+from a type that stores the problem — and a wrapper must store it, because
+`solve(&mut self) -> Solution<'_, ..>` hands back a solution that borrows from it,
+so the problem has to live somewhere that outlives the solution, and a field is
+the only such place. Such a wrapper could not then register a callback of its own:
+
+```rust
+struct Wrapper<'a, ..> {
+    my_callback: MyCallback,           // owned here
+    problem: Option<Problem<'a, ..>>,  // and borrowed from, here
+}
+```
+
+That is self-referential, and moving the wrapper would dangle the reference, so
+the compiler is right to reject it. The workaround in the dependent crate was a
+zero-sized callback, leaked to get a `&'static mut`.
+
+`with_owned_callback(Box<dyn Callback<T, D, P> + 'o>)` is the addition. Storage is
+one list of either kind:
+
+```rust
+enum Registered<'a, 'o, T: Scalar, const D: usize, P> {
+    Borrowed(&'a mut dyn Callback<T, D, P>),
+    Owned(Box<dyn Callback<T, D, P> + 'o>),
+}
+```
+
+`Registered` implements `Callback` itself, forwarding both methods, which is why
+none of the event handling in `solve` changed: the immutable scan for the earliest
+crossing goes through `event` and the single effect that fires goes through
+`effect`, whichever kind is holding it. Root finding, jump recording, the
+controller reset and `EVENT_PROBE_FRACTION` are untouched.
+
+### The second lifetime, which was not in the plan
+
+`Problem` gained a lifetime parameter: `Problem<'a, 'o, T, D, S, P, C>`. This was
+meant to be purely additive and it cannot be, for a reason worth writing down.
+
+Dropping a `Box<dyn Callback>` runs a destructor the compiler cannot see into, so
+dropck requires everything that box's lifetime covers to be alive when the problem
+is dropped. With one lifetime for both kinds, that requirement lands on the
+borrowed callbacks too — and it is type-based, so it applies even to a problem that
+owns no callback at all. The first attempt, a single `'a`, failed exactly here:
+
+```
+error[E0502]: cannot borrow `once.fired` as immutable because it is also borrowed as mutable
+  |
+  | mutable borrow might be used here, when `problem` is dropped
+```
+
+That is `test_callback_can_remember_firing`, reading a borrowed callback's own
+state after the solve, which is the only way a caller learns what a stateful
+callback did and the whole reason the borrowed form exists. Keeping it means the
+owned box needs a lifetime of its own, so that only `'o` has to outlive the drop
+and NLL can still end a borrowed callback's borrow at the problem's last use.
+
+The cost is arity: code that names the type explicitly needs
+`Problem<'a, '_, ..>`. Both lifetimes are inferred, so nothing that only calls
+`Problem::new(..).with_callback(..).solve()` is affected.
+
+### Not done: getting owned callbacks back
+
+Deliberately skipped. Downcasting a `Box<dyn Callback>` to its concrete type needs
+`Any`, which means a `'static` bound and an `as_any` method on every implementor —
+a cost paid by all callbacks for a need that has two cheaper answers already: stay
+borrowed, or share an `Rc<Cell<_>>`, which can be read while the solution is still
+alive.
+
+### Measurement
+
+The benches register no callbacks, so what they exercise is the list's element type
+changing from a reference to a two-variant enum, on a list that is empty. Measured
+by alternating old and new rather than running one after the other, because the
+first attempt at it reported a 4% improvement, which nothing in this change could
+produce: between-sitting drift on this machine is larger than the effect being
+looked for.
+
+| bench | old | new |
+|---|---|---|
+| `bench_orbit` | 356.06, 358.91 us | 357.98, 355.18 us |
+| `bench_simple_1d` | 877.7, 881.1 ns | 877.4, 895.7 ns |
+| `bench_interpolation_1d` | 2.2683, 2.2104 us | 2.2444, 2.2496 us |
+
+Two runs of the same binary differ by up to 2.6% here, and the old and new ranges
+overlap on all three. No measurable change.
+
+### Testing
+
+- an owned callback producing the same times and states as the identical borrowed
+  one, element for element;
+- an owned stateful callback escalating on its second firing, read off the
+  trajectory since its state cannot be read directly;
+- owned and borrowed callbacks on one problem with their events at the same time,
+  where the final state distinguishes registration order, and the borrowed one's
+  count is read after the solve;
+- a boxed `FnCallback`, which holds two closure references, confirming the `+ 'o`
+  rather than `'static`;
+- the wrapper from the top of this section, compiling and solving.
+
 ## Known Limitations
 
 - An even number of crossings inside one step is still missed: the two ends of

src/controller.rs

@@ -76,9 +76,7 @@ impl<T: Scalar> Controller<T> for PIController<T> {
             TryStep::Accepted(prev_step, h)
         } else {
             // Reject the stepsize and propose a smaller one for the current step
-            TryStep::NotYetAccepted(
-                prev_step / self.factor_c1.min(factor_11 / self.safety_factor),
-            )
+            TryStep::NotYetAccepted(prev_step / self.factor_c1.min(factor_11 / self.safety_factor))
         }
     }
 
@@ -124,7 +122,6 @@ impl<T: Scalar> Default for PIController<T> {
     }
 }
 
-
 /// PID (Proportional-Integral-Derivative) step size controller
 ///
 /// The PID controller is an advanced adaptive time-stepping controller that provides
@@ -171,20 +168,20 @@ impl<T: Scalar> Default for PIController<T> {
 #[derive(Debug, Clone, Copy)]
 pub struct PIDController<T: Scalar = f64> {
     // PID Coefficients
-    pub p_term: f64,  // Proportional: reaction to current error
-    pub i_term: f64,  // Integral: reaction to error history
-    pub d_term: f64,  // Derivative: reaction to error rate of change
+    pub p_term: f64, // Proportional: reaction to current error
+    pub i_term: f64, // Integral: reaction to error history
+    pub d_term: f64, // Derivative: reaction to error rate of change
 
     // Constraints
-    pub factor_c1: f64,      // 1/min_factor (maximum step decrease)
-    pub factor_c2: f64,      // 1/max_factor (maximum step increase)
-    pub h_max: f64,          // Maximum allowed step size
-    pub safety_factor: f64,  // Safety factor (typically 0.9)
+    pub factor_c1: f64,     // 1/min_factor (maximum step decrease)
+    pub factor_c2: f64,     // 1/max_factor (maximum step increase)
+    pub h_max: f64,         // Maximum allowed step size
+    pub safety_factor: f64, // Safety factor (typically 0.9)
 
     // Error history for PID control
-    pub err_old: f64,     // ε_{n-1}: previous step error
-    pub err_older: f64,   // ε_{n-2}: error two steps ago
-    pub h_old: f64,       // h_{n-1}: previous step size
+    pub err_old: f64,   // ε_{n-1}: previous step error
+    pub err_older: f64, // ε_{n-2}: error two steps ago
+    pub h_old: f64,     // h_{n-1}: previous step size
 
     // Next step guess
     pub next_step_guess: TryStep<T>,
@@ -199,15 +196,15 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
         let factor = if self.err_old <= 0.0 {
             // First step: use only proportional control
             let factor_11 = err.powf(self.p_term);
-            self.factor_c2.max(
-                self.factor_c1.min(factor_11 / self.safety_factor)
-            )
+            self.factor_c2
+                .max(self.factor_c1.min(factor_11 / self.safety_factor))
         } else if self.err_older <= 0.0 {
             // Second step: use PI control (proportional + integral)
             let factor_11 = err.powf(self.p_term);
             let factor_12 = self.err_old.powf(-self.i_term);
             self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 / self.safety_factor)
+                self.factor_c1
+                    .min(factor_11 * factor_12 / self.safety_factor),
             )
         } else {
             // Full PID control (proportional + integral + derivative)
@@ -220,7 +217,8 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
                 1.0
             };
             self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 * factor_13 / self.safety_factor)
+                self.factor_c1
+                    .min(factor_11 * factor_12 * factor_13 / self.safety_factor),
             )
         };
 
@@ -230,7 +228,7 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
 
             // Update error history for next step
             self.err_older = self.err_old;
-            self.err_old = err.max(1.0e-4);  // Prevent very small values
+            self.err_old = err.max(1.0e-4); // Prevent very small values
             self.h_old = h.re();
 
             // Apply maximum step size limit
@@ -297,9 +295,9 @@ impl<T: Scalar> PIDController<T> {
             factor_c2: 1.0 / max_factor,
             h_max: h_max.abs(),
             safety_factor,
-            err_old: 0.0,      // No history initially
-            err_older: 0.0,    // No history initially
-            h_old: 0.0,        // No history initially
+            err_old: 0.0,   // No history initially
+            err_older: 0.0, // No history initially
+            h_old: 0.0,     // No history initially
             next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)),
         }
     }
@@ -317,14 +315,14 @@ impl<T: Scalar> PIDController<T> {
     pub fn for_order(order: usize) -> Self {
         let k_plus_1 = (order + 1) as f64;
         Self::new(
-            0.49 / k_plus_1,  // p_term: proportional
-            0.34 / k_plus_1,  // i_term: integral
-            0.10 / k_plus_1,  // d_term: derivative
-            10.0,              // max_factor
-            0.2,               // min_factor
-            100000.0,          // h_max
-            0.9,               // safety_factor
-            1e-4,              // initial_h
+            0.49 / k_plus_1, // p_term: proportional
+            0.34 / k_plus_1, // i_term: integral
+            0.10 / k_plus_1, // d_term: derivative
+            10.0,            // max_factor
+            0.2,             // min_factor
+            100000.0,        // h_max
+            0.9,             // safety_factor
+            1e-4,            // initial_h
         )
     }
 }
@@ -342,19 +340,18 @@ impl<T: Scalar> Default for PIDController<T> {
     /// For method-specific tuning, use `PIDController::for_order(order)` instead.
     fn default() -> Self {
         Self::new(
-            0.07,      // p_term (proportional)
-            0.04,      // i_term (integral)
-            0.01,      // d_term (derivative)
-            10.0,      // max_factor
-            0.2,       // min_factor
-            100000.0,  // h_max
-            0.9,       // safety_factor
-            1e-4,      // initial_h
+            0.07,     // p_term (proportional)
+            0.04,     // i_term (integral)
+            0.01,     // d_term (derivative)
+            10.0,     // max_factor
+            0.2,      // min_factor
+            100000.0, // h_max
+            0.9,      // safety_factor
+            1e-4,     // initial_h
         )
     }
 }
 
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -406,7 +403,7 @@ mod tests {
         // Default uses conservative PID coefficients
         assert_eq!(controller.p_term, 0.07);
         assert_eq!(controller.i_term, 0.04);
-        assert_eq!(controller.d_term, 0.01);  // True PID with derivative term
+        assert_eq!(controller.d_term, 0.01); // True PID with derivative term
     }
 
     #[test]
@@ -443,7 +440,8 @@ mod tests {
 
         for &err in &errors {
             let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
+            let mut pid_result =
+                <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
 
             if pi_result.is_accepted() {
                 pi_steps.push(pi_result.extract());
@@ -555,7 +553,8 @@ mod tests {
 
         for &err in &errors {
             let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
+            let mut pid_result =
+                <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
 
             if pi_result.is_accepted() {
                 let pi_h_next = pi_result.reset().unwrap().extract();
@@ -574,11 +573,8 @@ mod tests {
 
         // Compute standard deviation of log step size ratios
         let pi_mean: f64 = pi_ratios.iter().sum::<f64>() / pi_ratios.len() as f64;
-        let pi_variance: f64 = pi_ratios
-            .iter()
-            .map(|r| (r - pi_mean).powi(2))
-            .sum::<f64>()
-            / pi_ratios.len() as f64;
+        let pi_variance: f64 =
+            pi_ratios.iter().map(|r| (r - pi_mean).powi(2)).sum::<f64>() / pi_ratios.len() as f64;
         let pi_std = pi_variance.sqrt();
 
         let pid_mean: f64 = pid_ratios.iter().sum::<f64>() / pid_ratios.len() as f64;

src/integrator/bs3.rs

@@ -64,7 +64,7 @@ pub struct BS3<const D: usize> {
 impl<const D: usize> Default for BS3<D>
 where
     BS3<D>: Integrator<D>,
- {
+{
     fn default() -> Self {
         Self::new()
     }
@@ -120,6 +120,10 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     // Row 2: [1/2]
     // Row 3: [0, 3/4]
     // Row 4: [2/9, 1/3, 4/9]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
     const A: &'a [f64] = &[
         1.0 / 2.0,           // a[1,0]
         0.0,                 // a[2,0]
@@ -130,6 +134,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Solution weights (3rd order)
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const B: &'a [f64] = &[
         2.0 / 9.0,           // b[0]
         1.0 / 3.0,           // b[1]
@@ -138,6 +144,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Error estimate weights (difference between 3rd and 2nd order)
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const B_ERROR: &'a [f64] = &[
         2.0 / 9.0 - 7.0 / 24.0,      // b[0] - b*[0]
         1.0 / 3.0 - 1.0 / 4.0,       // b[1] - b*[1]
@@ -146,6 +154,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Stage times
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const C: &'a [f64] = &[
         0.0,                 // c[0]
         1.0 / 2.0,           // c[1]
@@ -191,7 +201,9 @@ where
 
         // Stage 4: k4 = f(t + c[3]*h, y + h*(a[3,0]*k1 + a[3,1]*k2 + a[3,2]*k3))
         let y4 = ode.y
-            + (k[0] * T::lift(Self::A[3]) + k[1] * T::lift(Self::A[4]) + k[2] * T::lift(Self::A[5]))
+            + (k[0] * T::lift(Self::A[3])
+                + k[1] * T::lift(Self::A[4])
+                + k[2] * T::lift(Self::A[5]))
                 * h;
         k[3] = (ode.f)(ode.t + h * Self::C[3], y4, &ode.params);
 
@@ -225,10 +237,10 @@ where
         // BS3 uses standard cubic Hermite interpolation with derivatives at endpoints
         // Store: y0, y1, f0=k[0], f1=k[3] (FSAL)
         let dense_coeffs = vec![
-            ode.y,           // y0 at start of step
-            next_y,          // y1 at end of step
-            k[0],            // f(t0, y0) - derivative at start
-            k[3],            // f(t1, y1) - derivative at end (FSAL)
+            ode.y,  // y0 at start of step
+            next_y, // y1 at end of step
+            k[0],   // f(t0, y0) - derivative at start
+            k[3],   // f(t1, y1) - derivative at end (FSAL)
         ];
 
         StepResult {
@@ -304,7 +316,7 @@ mod tests {
         let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
 
         let bs3 = BS3::new();
-        let h = 0.001;  // Smaller step size for tighter tolerances
+        let h = 0.001; // Smaller step size for tighter tolerances
 
         let __step = bs3.step(&ode, h);
         let (y_next, err, dense) = (__step.y, __step.err, __step.dense);
@@ -334,7 +346,7 @@ mod tests {
         let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
 
         let bs3 = BS3::new();
-        let h = 0.001;  // Smaller step size
+        let h = 0.001; // Smaller step size
 
         let __step = bs3.step(&ode, h);
         let (_y_next, _err, dense) = (__step.y, __step.err, __step.dense);

src/integrator/dormand_prince.rs

@@ -212,7 +212,11 @@ where
             for (j, item) in k.iter().enumerate().take(i) {
                 y_term += item * T::lift(Self::A[(i * (i - 1)) / 2 + j]);
             }
-            k[i] = (ode.f)(ode.t + T::lift(Self::C[i]) * h, ode.y + y_term * h, &ode.params);
+            k[i] = (ode.f)(
+                ode.t + T::lift(Self::C[i]) * h,
+                ode.y + y_term * h,
+                &ode.params,
+            );
 
             // Use that and bis to calculate the y and error terms
             next_y += k[i] * h * T::lift(Self::B[i]);

src/integrator/method.rs

@@ -406,7 +406,6 @@ mod tests {
     /// integrator it wraps does, which is only true if nothing erases the
     /// scalar. Vern7 is the demanding case: interpolating inside a step calls
     /// the right-hand side again, in dual arithmetic, after the solve.
-    #[cfg(feature = "autodiff")]
     #[test]
     fn test_method_preserves_autodiff() {
         use num_dual::Dual64;

src/integrator/rosenbrock.rs

@@ -28,10 +28,11 @@ impl Default for JacobianUpdateStrategy {
 
 /// Compute the Jacobian matrix ∂f/∂y using forward finite differences
 ///
-/// For a system y' = f(t, y), computes the D×D Jacobian matrix J where J[i,j] = ∂f_i/∂y_j
+/// For a system `y' = f(t, y)`, computes the D×D Jacobian matrix `J` where
+/// `J[i,j] = ∂f_i/∂y_j`
 ///
-/// Uses forward differences: J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε
-/// where ε = √(machine_epsilon) * max(|y[j]|, 1.0)
+/// Uses forward differences: `J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε`
+/// where `ε = √(machine_epsilon) * max(|y[j]|, 1.0)`
 pub fn compute_jacobian<T: Scalar, const D: usize, P>(
     f: ProblemFunction<'_, T, D, P>,
     t: T,
@@ -443,12 +444,21 @@ mod tests {
         // So log(error) = 2*log(h) + const
         // Slope should be approximately 2
         for i in 0..errors.len() - 1 {
-            let rate =
-                (errors[i].log10() - errors[i + 1].log10()) / (step_sizes[i].log10() - step_sizes[i + 1].log10());
-            println!("Convergence rate between h={} and h={}: {}", step_sizes[i], step_sizes[i+1], rate);
+            let rate = (errors[i].log10() - errors[i + 1].log10())
+                / (step_sizes[i].log10() - step_sizes[i + 1].log10());
+            println!(
+                "Convergence rate between h={} and h={}: {}",
+                step_sizes[i],
+                step_sizes[i + 1],
+                rate
+            );
 
             // Should be close to 2 for a 2nd order method
-            assert!(rate > 1.8 && rate < 2.2, "Convergence rate {} is not close to 2", rate);
+            assert!(
+                rate > 1.8 && rate < 2.2,
+                "Convergence rate {} is not close to 2",
+                rate
+            );
         }
     }
 
@@ -466,7 +476,7 @@ mod tests {
         let analytical = |t: f64| 0.5 * (1.01 * t).exp();
 
         // Test convergence with Julia's step sizes: (1/2)^(6:-1:3) = [1/64, 1/32, 1/16, 1/8]
-        let step_sizes = vec![1.0/64.0, 1.0/32.0, 1.0/16.0, 1.0/8.0];
+        let step_sizes = vec![1.0 / 64.0, 1.0 / 32.0, 1.0 / 16.0, 1.0 / 8.0];
         let mut errors = Vec::new();
 
         for &h in &step_sizes {
@@ -504,8 +514,11 @@ mod tests {
         println!("Orders per step refinement: {:?}", orders);
 
         // Julia tests: @test sim.𝒪est[:final]≈2 atol=0.2
-        assert!((avg_order - 2.0).abs() < 0.2,
-                "Convergence order {:.2} not within 0.2 of expected order 2", avg_order);
+        assert!(
+            (avg_order - 2.0).abs() < 0.2,
+            "Convergence order {:.2} not within 0.2 of expected order 2",
+            avg_order
+        );
     }
 
     #[test]
@@ -528,11 +541,16 @@ mod tests {
         let mut problem = Problem::new(ode, rb23, controller);
         let solution = problem.solve();
 
-        println!("Adaptive solve completed in {} steps", solution.states.len());
+        println!(
+            "Adaptive solve completed in {} steps",
+            solution.states.len()
+        );
 
-        assert!(solution.states.len() < 20,
-                "Adaptive solve should complete in less than 20 steps, got {}",
-                solution.states.len());
+        assert!(
+            solution.states.len() < 20,
+            "Adaptive solve should complete in less than 20 steps, got {}",
+            solution.states.len()
+        );
 
         // Verify final value is accurate
         let analytical = 0.5 * (1.01_f64 * 1.0).exp();

src/integrator/vern7.rs

@@ -80,7 +80,7 @@ pub struct Vern7<const D: usize> {
 impl<const D: usize> Default for Vern7<D>
 where
     Vern7<D>: Integrator<D>,
- {
+{
     fn default() -> Self {
         Self::new()
     }
@@ -141,6 +141,8 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     // Stage 8: [a81, 0, a83, a84, a85, a86, a87]
     // Stage 9: [a91, 0, a93, a94, a95, a96, a97, a98]
     // Stage 10: [a101, 0, a103, a104, a105, a106, a107, 0, 0]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
     const A: &'a [f64] = &[
         // Stage 2
         0.005,
@@ -177,6 +179,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     ];
 
     // Error estimate weights (difference between 7th and 6th order: b - b*)
+    #[rustfmt::skip]
     const B_ERROR: &'a [f64] = &[
         0.002547011879931045,   // b1 - b*1
         0.0,                    // b2 - b*2
@@ -191,6 +194,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     ];
 
     // Time nodes (c coefficients)
+    #[rustfmt::skip]
     const C: &'a [f64] = &[
         0.0,                  // c1
         0.005,                // c2
@@ -210,6 +214,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
 
 impl<'a, const D: usize> Vern7ExtraStages<'a> for Vern7<D> {
     // Time nodes for extra stages
+    #[rustfmt::skip]
     const C_EXTRA: &'a [f64] = &[
         1.0,    // c11
         0.29,   // c12
@@ -375,6 +380,8 @@ where
         }
     }
 
+    // One interpolation coefficient per line, likewise.
+    #[rustfmt::skip]
     fn interpolate<T: Scalar>(
         &self,
         t_start: T,
@@ -744,7 +751,7 @@ mod tests {
         }
 
         let y0 = Vector1::new(1.0);
-        let t_end: f64 = 1.0;  // Longer interval to get larger errors
+        let t_end: f64 = 1.0; // Longer interval to get larger errors
         let exact = t_end.exp();
 
         let step_sizes: [f64; 3] = [0.2, 0.1, 0.05];
@@ -771,16 +778,8 @@ mod tests {
         let ratio2 = errors[1] / errors[2];
 
         // Allow some tolerance (expect ratio between 64 and 256)
-        assert!(
-            ratio1 > 64.0 && ratio1 < 256.0,
-            "First ratio: {}",
-            ratio1
-        );
-        assert!(
-            ratio2 > 64.0 && ratio2 < 256.0,
-            "Second ratio: {}",
-            ratio2
-        );
+        assert!(ratio1 > 64.0 && ratio1 < 256.0, "First ratio: {}", ratio1);
+        assert!(ratio2 > 64.0 && ratio2 < 256.0, "Second ratio: {}", ratio2);
     }
 
     #[test]
@@ -828,7 +827,7 @@ mod tests {
             Vector2::new(y[1], -y[0])
         }
 
-        let y0 = Vector2::new(1.0, 0.0);  // Start at maximum displacement, zero velocity
+        let y0 = Vector2::new(1.0, 0.0); // Start at maximum displacement, zero velocity
 
         // Period of harmonic oscillator is 2π
         let period = 2.0 * std::f64::consts::PI;
@@ -854,7 +853,10 @@ mod tests {
 
         println!("Initial energy: {}", energy_initial);
         println!("Final energy: {}", energy_final);
-        println!("Energy drift after {} periods: {:.2e}", num_periods, energy_drift);
+        println!(
+            "Energy drift after {} periods: {:.2e}",
+            num_periods, energy_drift
+        );
         println!("Number of steps: {}", solution.times.len());
 
         // Energy should be conserved to high precision (< 1e-7 relative error over 1000 periods)

src/lib.rs

@@ -1,8 +1,7 @@
 #![allow(dead_code)]
 // The readme is the crate's front page, and including it here means every
-// example in it is compiled and run by `cargo test`. Requires `autodiff`
-// because one of those examples differentiates.
-#![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))]
+// example in it is compiled and run by `cargo test`.
+#![doc = include_str!("../readme.md")]
 
 pub mod callback;
 pub mod controller;
@@ -19,8 +18,8 @@ pub mod prelude {
     pub use super::integrator::method::{Method, MethodName};
     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::ode::ODE;
     pub use super::problem::{Problem, Solution};
     pub use super::scalar::Scalar;
 }
@@ -160,7 +159,7 @@ mod tests {
 }
 
 /// Tests that the derivatives a dual instantiation produces are the right ones.
-#[cfg(all(test, feature = "autodiff"))]
+#[cfg(test)]
 mod autodiff_tests {
     use crate::ode::ProblemFunction;
     use crate::prelude::*;
@@ -376,7 +375,13 @@ mod autodiff_tests {
         }
 
         let y0 = Dual64::from_re(1.0).derivative();
-        let ode = ODE::new(&derivative, Dual64::lift(0.0), Dual64::lift(2.0), Vector1::new(y0), ());
+        let ode = ODE::new(
+            &derivative,
+            Dual64::lift(0.0),
+            Dual64::lift(2.0),
+            Vector1::new(y0),
+            (),
+        );
         let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
         let mut problem = Problem::new(ode, vern7, PIController::default());
         let solution = problem.solve();

src/problem.rs

@@ -33,7 +33,56 @@ use super::scalar::Scalar;
 /// the event.
 const EVENT_PROBE_FRACTION: f64 = 1e-6;
 
-pub struct Problem<'a, T: Scalar, const D: usize, S, P, C>
+/// A callback a [`Problem`] watches, which the problem either borrows or owns.
+///
+/// Borrowing is what lets a caller read a stateful callback's own fields after
+/// the solve; owning is what lets a struct that stores the problem register a
+/// callback of its own, which it cannot do by reference without becoming
+/// self-referential.
+///
+/// Implementing [`Callback`] on this, rather than reaching inside it, is what
+/// keeps the event handling in [`Problem::solve`] written once: the immutable
+/// scan for a crossing goes through `event` and the one effect that fires goes
+/// through `effect`, whichever kind is holding it.
+///
+/// The two lifetimes are separate on purpose. Dropping a `Box<dyn Callback>` runs
+/// a destructor the compiler cannot see into, so everything that box's lifetime
+/// covers has to still be alive when the problem is dropped. Were that the same
+/// lifetime as the borrowed callbacks', a borrowed callback could not be read
+/// back until the problem went out of scope, and reading one back is the whole
+/// reason the borrowed form exists.
+enum Registered<'a, 'o, T: Scalar, const D: usize, P> {
+    Borrowed(&'a mut dyn Callback<T, D, P>),
+    Owned(Box<dyn Callback<T, D, P> + 'o>),
+}
+
+impl<T: Scalar, const D: usize, P> Callback<T, D, P> for Registered<'_, '_, T, D, P> {
+    #[inline]
+    fn event(&self, t: T, y: &SVector<T, D>, params: &P) -> T {
+        match self {
+            Registered::Borrowed(callback) => callback.event(t, y, params),
+            Registered::Owned(callback) => callback.event(t, y, params),
+        }
+    }
+
+    #[inline]
+    fn effect(&mut self, ode: &mut ODE<'_, T, D, P>) {
+        match self {
+            Registered::Borrowed(callback) => callback.effect(ode),
+            Registered::Owned(callback) => callback.effect(ode),
+        }
+    }
+}
+
+/// An ODE, an integrator, a controller and any callbacks, ready to solve.
+///
+/// The first lifetime is the problem's own borrows: the right-hand side, and any
+/// callback registered with [`with_callback`](Problem::with_callback). The second
+/// is what an owned callback from
+/// [`with_owned_callback`](Problem::with_owned_callback) is allowed to borrow.
+/// Both are inferred, so `Problem<'_, '_, f64, 3, _, _, _>` is the usual way to
+/// name the type when it has to be named at all.
+pub struct Problem<'a, 'o, T: Scalar, const D: usize, S, P, C>
 where
     S: Integrator<D>,
     C: Controller<T>,
@@ -41,10 +90,10 @@ where
     ode: ODE<'a, T, D, P>,
     integrator: S,
     controller: C,
-    callbacks: Vec<&'a mut dyn Callback<T, D, P>>,
+    callbacks: Vec<Registered<'a, 'o, T, D, P>>,
 }
 
-impl<'a, T: Scalar, const D: usize, S, P, C> Problem<'a, T, D, S, P, C>
+impl<'a, 'o, T: Scalar, const D: usize, S, P, C> Problem<'a, 'o, T, D, S, P, C>
 where
     S: Integrator<D> + Copy,
     C: Controller<T>,
@@ -91,8 +140,7 @@ where
         // crossing starts from rather than at — see `EVENT_PROBE_FRACTION`.
         let mut last_fired: Vec<Option<f64>> = vec![None; self.callbacks.len()];
         while !self.ode.finished && self.ode.t.re() < self.ode.t_end.re() {
-            if (self.ode.t + self.controller.next_step_guess().extract()).re()
-                > self.ode.t_end.re()
+            if (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
@@ -149,7 +197,9 @@ where
                 // The event condition is a real-valued question, so it is
                 // asked in the reals — see this module's header.
                 let before = if search_start.re() == 0.0 {
-                    callback.event(self.ode.t, &self.ode.y, &self.ode.params).re()
+                    callback
+                        .event(self.ode.t, &self.ode.y, &self.ode.params)
+                        .re()
                 } else {
                     let probe = self.integrator.step(&self.ode, search_start).y;
                     callback
@@ -171,7 +221,9 @@ where
                     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()
+                        callback
+                            .event(self.ode.t + h, &test_y, &self.ode.params)
+                            .re()
                     };
                     Some(
                         find_root_brent(search_start.re(), curr_step.re(), &f, &mut convergency)
@@ -235,12 +287,88 @@ where
         }
     }
 
-    /// Registers a callback.
+    /// Registers a borrowed callback.
     ///
     /// Callbacks are watched in the order they are registered, which decides
-    /// which one fires when two events land at the same time.
+    /// which one fires when two events land at the same time, whether they are
+    /// borrowed or owned.
+    ///
+    /// Borrowing means something else owns the callback and outlives the problem,
+    /// and it is what lets a caller read the callback's own state once the solve
+    /// is done:
+    ///
+    /// ```
+    /// # use nalgebra::Vector1;
+    /// # use ordinary_diffeq::prelude::*;
+    /// struct Count {
+    ///     fired: usize,
+    /// }
+    ///
+    /// impl Callback<f64, 1, ()> for Count {
+    ///     fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
+    ///         1.0 - t
+    ///     }
+    ///
+    ///     fn effect(&mut self, _ode: &mut ODE<'_, f64, 1, ()>) {
+    ///         self.fired += 1;
+    ///     }
+    /// }
+    ///
+    /// fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+    ///     Vector1::new(1.0)
+    /// }
+    ///
+    /// let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(0.0), ());
+    /// let mut counter = Count { fired: 0 };
+    /// let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+    ///     .with_callback(&mut counter);
+    /// problem.solve();
+    ///
+    /// // The borrow ended with the problem's last use, so the state is readable.
+    /// assert_eq!(counter.fired, 1);
+    /// ```
     pub fn with_callback(mut self, callback: &'a mut dyn Callback<T, D, P>) -> Self {
-        self.callbacks.push(callback);
+        self.callbacks.push(Registered::Borrowed(callback));
+        self
+    }
+
+    /// Registers a callback the problem owns.
+    ///
+    /// For a caller that cannot lend one: a struct that stores the problem in a
+    /// field — which it must, since [`solve`](Problem::solve) returns a solution
+    /// borrowing from it — cannot also hold the callback in a field and lend it
+    /// to the problem, because that makes the struct self-referential and moving
+    /// it would dangle the reference. Moving the callback into the problem avoids
+    /// the question.
+    ///
+    /// The `+ 'a` is deliberate: an owned callback may hold borrows of its own,
+    /// which is exactly what [`FnCallback`](crate::callback::FnCallback) is.
+    ///
+    /// Nothing can read an owned callback's state afterwards. A callback that has
+    /// something to report either stays borrowed, or shares a handle — an
+    /// `Rc<Cell<_>>` — with its owner, which works while the solution is still
+    /// alive.
+    ///
+    /// ```
+    /// # use nalgebra::Vector1;
+    /// # use ordinary_diffeq::prelude::*;
+    /// # fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
+    /// #     Vector1::new(1.0)
+    /// # }
+    /// let ode = ODE::new(&derivative, 0.0, 10.0, Vector1::new(0.0), ());
+    /// let jump = FnCallback {
+    ///     event: &|_t: f64, y: &Vector1<f64>, _p: &()| 5.0 - y[0],
+    ///     effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
+    /// };
+    ///
+    /// let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+    ///     .with_owned_callback(Box::new(jump));
+    /// let solution = problem.solve();
+    ///
+    /// assert!((solution.states.last().unwrap()[0] - 105.0).abs() < 1e-9);
+    /// ```
+    pub fn with_owned_callback(mut self, callback: Box<dyn Callback<T, D, P> + 'o>) -> Self {
+        self.callbacks.push(Registered::Owned(callback));
         self
     }
 }
@@ -366,6 +494,7 @@ mod tests {
     use crate::callback::{stop, FnCallback};
     use crate::controller::PIController;
     use crate::integrator::dormand_prince::DormandPrince45;
+    use crate::integrator::vern7::Vern7;
     use crate::integrator::Tolerances;
     use approx::assert_relative_eq;
     use nalgebra::{Matrix3, Vector1, Vector2, Vector3};
@@ -440,8 +569,7 @@ mod tests {
             effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
         };
 
-        let mut problem =
-            Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
         let solution = problem.solve();
 
         // Five units of integration before the jump, five after.
@@ -464,8 +592,7 @@ mod tests {
             effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
         };
 
-        let mut problem =
-            Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
         let solution = problem.solve();
 
         let at_crossing: Vec<usize> = (0..solution.times.len())
@@ -650,8 +777,7 @@ mod tests {
             effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = -ode.y,
         };
 
-        let mut problem =
-            Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
         let solution = problem.solve();
 
         // y rises to 1 at t = 1 and flips to -1, which the effect leaves on the
@@ -665,7 +791,11 @@ mod tests {
         assert_relative_eq!(jumps[0], 1.0, max_relative = 1e-9);
         assert_relative_eq!(jumps[1], 3.0, max_relative = 1e-9);
         assert_relative_eq!(solution.times.last().unwrap(), &3.0, max_relative = 1e-12);
-        assert_relative_eq!(solution.states.last().unwrap()[0], -1.0, max_relative = 1e-8);
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            -1.0,
+            max_relative = 1e-8
+        );
     }
 
     /// The controller's error history is discarded when an effect changes the
@@ -708,8 +838,8 @@ mod tests {
                 effect,
             };
             let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
-            let mut problem = Problem::new(ode, DormandPrince45::new(), spy)
-                .with_callback(&mut callback);
+            let mut problem =
+                Problem::new(ode, DormandPrince45::new(), spy).with_callback(&mut callback);
             problem.solve();
             resets.get()
         }
@@ -722,6 +852,187 @@ mod tests {
         assert_eq!(resets_during(&|_ode: &mut ODE<'_, f64, 1, ()>| {}), 0);
     }
 
+    /// An owned callback does what the same callback does when borrowed: same
+    /// crossings, same jumps, same trajectory.
+    #[test]
+    fn test_owned_callback_matches_borrowed() {
+        fn jump_at_five<'c>() -> FnCallback<'c, f64, 1, ()> {
+            FnCallback {
+                event: &|_t: f64, y: &Vector1<f64>, _p: &()| 5.0 - y[0],
+                effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
+            }
+        }
+
+        let borrowed = {
+            let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+            let mut callback = jump_at_five();
+            let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+                .with_callback(&mut callback);
+            let solution = problem.solve();
+            (solution.times.clone(), solution.states.clone())
+        };
+
+        let owned = {
+            let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+            let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+                .with_owned_callback(Box::new(jump_at_five()));
+            let solution = problem.solve();
+            (solution.times.clone(), solution.states.clone())
+        };
+
+        assert_eq!(owned.0, borrowed.0);
+        assert_eq!(owned.1, borrowed.1);
+        assert_relative_eq!(owned.1.last().unwrap()[0], 105.0, max_relative = 1e-9);
+    }
+
+    /// An owned callback keeps state across its firings, even though nothing can
+    /// read that state afterwards: what it did shows up in the trajectory.
+    #[test]
+    fn test_owned_stateful_callback() {
+        struct Escalate {
+            fired: usize,
+        }
+
+        impl Callback<f64, 1, ()> for Escalate {
+            fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
+                // Crosses at t = 1 and t = 3.
+                (t * PI / 2.0).cos()
+            }
+
+            fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) {
+                self.fired += 1;
+                ode.y += Vector1::new(if self.fired == 1 { 10.0 } else { 100.0 });
+            }
+        }
+
+        let ode = ODE::new(&unit_rate, 0.0, 4.0, Vector1::new(0.0), ());
+        // Capped so one step cannot swallow two crossings; see the borrowed
+        // equivalent in `test_callback_can_remember_firing`.
+        let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 1.0, 0.9, 1e-4);
+        let mut problem = Problem::new(ode, DormandPrince45::new(), controller)
+            .with_owned_callback(Box::new(Escalate { fired: 0 }));
+        let solution = problem.solve();
+
+        // 1 + 10, carried to 13, then + 100, carried to 114. Had the second
+        // firing added 10 as well it would end at 24.
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            114.0,
+            max_relative = 1e-9
+        );
+    }
+
+    /// Owned and borrowed callbacks on one problem. Both events land at the same
+    /// time, so registration order decides which effect applies first, and the
+    /// borrowed one's state is still readable once the solve is done.
+    #[test]
+    fn test_owned_and_borrowed_on_one_problem() {
+        struct AddOne {
+            fired: usize,
+        }
+
+        impl Callback<f64, 1, ()> for AddOne {
+            fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
+                2.0 - t
+            }
+
+            fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) {
+                self.fired += 1;
+                ode.y += Vector1::new(1.0);
+            }
+        }
+
+        let double = FnCallback {
+            event: &|t: f64, _y: &Vector1<f64>, _p: &()| 2.0 - t,
+            effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y *= 2.0,
+        };
+
+        let ode = ODE::new(&unit_rate, 0.0, 4.0, Vector1::new(0.0), ());
+        let mut add_one = AddOne { fired: 0 };
+        let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+            .with_owned_callback(Box::new(double))
+            .with_callback(&mut add_one);
+        let solution = problem.solve();
+        let last = solution.states.last().unwrap()[0];
+
+        // y = 2 at the crossing, doubled to 4 by the owned callback, then 5 from
+        // the borrowed one, then two more units of integration. The other order
+        // would give 8.
+        assert_relative_eq!(last, 7.0, max_relative = 1e-9);
+        assert_eq!(add_one.fired, 1);
+    }
+
+    /// An owned callback may hold borrows of its own, which is what `FnCallback`
+    /// is: a pair of closure references. The box is not `'static`.
+    #[test]
+    fn test_owned_callback_holding_a_borrow() {
+        let threshold = 5.0;
+        let target = Vector1::new(100.0);
+        let event = move |_t: f64, y: &Vector1<f64>, _p: &()| threshold - y[0];
+        let effect = move |ode: &mut ODE<'_, f64, 1, ()>| ode.y = target;
+
+        let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+        let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+            .with_owned_callback(Box::new(FnCallback {
+                event: &event,
+                effect: &effect,
+            }));
+        let solution = problem.solve();
+
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            105.0,
+            max_relative = 1e-9
+        );
+    }
+
+    /// The case owned callbacks exist for: a type that stores the problem,
+    /// because the solution borrows from it, and registers a callback of its own.
+    /// By reference that would make the wrapper self-referential and would not
+    /// compile.
+    #[test]
+    fn test_wrapper_owning_both_problem_and_callback() {
+        struct Shift;
+
+        impl Callback<f64, 1, ()> for Shift {
+            fn event(&self, _t: f64, y: &Vector1<f64>, _p: &()) -> f64 {
+                5.0 - y[0]
+            }
+
+            fn effect(&mut self, ode: &mut ODE<'_, f64, 1, ()>) {
+                ode.y = Vector1::new(100.0);
+            }
+        }
+
+        type Run<'a> = Problem<'a, 'a, f64, 1, DormandPrince45<1>, (), PIController<f64>>;
+
+        struct Wrapper<'a> {
+            problem: Option<Run<'a>>,
+        }
+
+        impl<'a> Wrapper<'a> {
+            fn solve(
+                &mut self,
+                ode: ODE<'a, f64, 1, ()>,
+            ) -> Solution<'_, f64, DormandPrince45<1>, 1, ()> {
+                let problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
+                    .with_owned_callback(Box::new(Shift));
+                self.problem = Some(problem);
+                self.problem.as_mut().unwrap().solve()
+            }
+        }
+
+        let mut wrapper = Wrapper { problem: None };
+        let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
+        let solution = wrapper.solve(ode);
+
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            105.0,
+            max_relative = 1e-9
+        );
+    }
+
     /// A bouncing ball, which is the standard exercise of repeated jumps: an
     /// event every time the ball reaches the floor, and an effect that reverses
     /// and shrinks its velocity.
@@ -739,7 +1050,7 @@ mod tests {
         }
 
         let ode = ODE::new(&free_fall, 0.0, 2.0, Vector2::new(1.0, 0.0), G);
-        let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
         let mut bounce = FnCallback {
             event: &|_t: f64, y: &Vector2<f64>, _p: &f64| y[0],
             effect: &|ode: &mut ODE<'_, f64, 2, f64>| {
@@ -803,7 +1114,7 @@ mod tests {
         let omega = Vector3::new(0.0, 0.0, 1.0);
         let t_end = 5.0; // past half a turn, short of a full one
         let ode = ODE::new(&mrp_rate, 0.0, t_end, Vector3::zeros(), omega);
-        let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
 
         let mut shadow = FnCallback {
             event: &|_t: f64, s: &Vector3<f64>, _p: &Vector3<f64>| 1.0 - s.norm(),
@@ -830,7 +1141,11 @@ mod tests {
         // And the attitude at the end is the one the body is actually in.
         let expected = -Vector3::new(0.0, 0.0, 1.0 / (t_end / 4.0).tan());
         for i in 0..3 {
-            assert_relative_eq!(solution.states.last().unwrap()[i], expected[i], epsilon = 1e-9);
+            assert_relative_eq!(
+                solution.states.last().unwrap()[i],
+                expected[i],
+                epsilon = 1e-9
+            );
         }
     }
 

src/scalar.rs

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