diff --git a/diffrax/_integrate.py b/diffrax/_integrate.py
index 5441e0a9..70d1abdc 100644
--- a/diffrax/_integrate.py
+++ b/diffrax/_integrate.py
@@ -455,9 +455,13 @@ def save_ts(subsaveat: SubSaveAt, save_state: SaveState) -> SaveState:
 
         def save_ts_impl(ts, fn, save_state: SaveState) -> SaveState:
             def _cond_fun(_save_state):
+                # A checkpoint recompute can round `tnext` two ULPs below the
+                # forward step end. A save that landed on that end still
+                # belongs to this step.
+                tnext = eqxi.nextafter(eqxi.nextafter(state.tnext))
                 return (
                     keep_step
-                    & (ts[_save_state.saveat_ts_index] <= state.tnext)
+                    & (ts[_save_state.saveat_ts_index] <= tnext)
                     & (_save_state.saveat_ts_index < len(ts))
                 )
 
diff --git a/test/test_adjoint.py b/test/test_adjoint.py
index 9e17e535..1c1b697a 100644
--- a/test/test_adjoint.py
+++ b/test/test_adjoint.py
@@ -453,3 +453,32 @@ def run_jvp(y0):
         return jax.jvp(run, (y0,), (jnp.ones_like(y0),))
 
     run_jvp(jnp.array(1.0))
+
+
+def test_recursive_checkpoint_adjoint_save_on_step_end():
+    # Save times that land on constant-step ends. See issue #788.
+    ts = jnp.linspace(0.0, 6.0, 11)
+    y0 = jnp.array([1.0, 0.0])
+
+    def loss(a, adjoint):
+        sol = diffrax.diffeqsolve(
+            diffrax.ODETerm(lambda t, y, a: jnp.stack([y[1], -a * y[0]])),
+            diffrax.Tsit5(),
+            0.0,
+            6.0,
+            0.6,
+            y0,
+            args=a,
+            saveat=diffrax.SaveAt(ts=ts),
+            stepsize_controller=diffrax.ConstantStepSize(),
+            adjoint=adjoint,
+            max_steps=20,
+        )
+        return jnp.sum(sol.ys[:, 0] * jnp.cos(ts))
+
+    a = jnp.array(0.7)
+    direct = jax.grad(lambda a: loss(a, diffrax.DirectAdjoint()))(a)
+    for checkpoints in (None, 4, 5, 6):
+        adjoint = diffrax.RecursiveCheckpointAdjoint(checkpoints=checkpoints)
+        got = jax.grad(lambda a: loss(a, adjoint))(a)
+        assert tree_allclose(got, direct)