Release 0.6.1

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

Cargo.toml

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

benches/bs3_vs_dp5.rs

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

benches/orbit.rs

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

benches/simple_1d.rs

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

benches/vern7_comparison.rs

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

hooks/pre-commit

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

hooks/pre-push

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

readme.md

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

roadmap/GETTING_STARTED.md

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

roadmap/README.md

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

roadmap/features/42-callback-effects.md

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

src/callback.rs

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

src/controller.rs

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

src/integrator/bs3.rs

@@ -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/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

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

src/ode.rs

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

src/problem.rs

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

src/scalar.rs

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