roadmap/README.md

@@ -298,6 +298,12 @@ Each feature below links to a detailed implementation plan in the `features/` di
   - **Dependencies**: None
   - **Effort**: Small
 
+- [ ] **[Reduce Per-Step Allocation](features/39-per-step-allocation.md)**
+  - Two heap allocations per accepted step (stage array, dense coefficients)
+  - Against ~150 ns of arithmetic; matters only for cheap right-hand sides
+  - **Dependencies**: None
+  - **Effort**: Medium
+
 ---
 
 ## Implementation Notes

roadmap/features/39-per-step-allocation.md

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

src/problem.rs

@@ -65,15 +65,10 @@ where
                 while !next_step_guess.is_accepted() {
                     // If that step isn't acceptable, then change the step until it is
                     result = self.integrator.step(&self.ode, next_step_guess.extract());
-                    // NOTE: `err` here is the error from the *previous* attempt, not
-                    // from the step just taken, and the assignment below only takes
-                    // effect on the next pass. This is preserved from 0.2.3 rather
-                    // than fixed, so that the 0.3.0 benchmarks measure the
-                    // genericisation and nothing else. See IMPROVEMENTS.md.
+                    err = result.err.unwrap();
                     next_step_guess = self
                         .controller
                         .determine_step(next_step_guess.extract(), err);
-                    err = result.err.unwrap();
                 }
                 // So at this point we can safely assume we have an accepted step
                 self.controller