diff --git a/diffrax/_solver/align.py b/diffrax/_solver/align.py
index c6bc6105..433b2779 100644
--- a/diffrax/_solver/align.py
+++ b/diffrax/_solver/align.py
@@ -6,6 +6,7 @@
 
 from .._custom_types import (
     AbstractSpaceTimeLevyArea,
+    Args,
     RealScalarLike,
 )
 from .._local_interpolation import LocalLinearInterpolation
@@ -156,6 +157,7 @@ def _compute_step(
         coeffs: _ALIGNCoeffs,
         rho: UnderdampedLangevinX,
         prev_f: UnderdampedLangevinX,
+        args: Args,
     ) -> tuple[
         UnderdampedLangevinX,
         UnderdampedLangevinX,
@@ -176,7 +178,7 @@ def _compute_step(
             - coeffs.b1**ω * uh**ω * f0**ω
             + rho**ω * (coeffs.b1**ω * w**ω + coeffs.chh**ω * hh**ω)
         ).ω
-        f1 = f(x1)
+        f1 = f(x1, args)
         v1 = (
             coeffs.beta**ω * v0**ω
             - u**ω * ((coeffs.a1**ω - coeffs.b1**ω) * f0**ω + coeffs.b1**ω * f1**ω)
diff --git a/diffrax/_solver/foster_langevin_srk.py b/diffrax/_solver/foster_langevin_srk.py
index dbdf3939..47ae3090 100644
--- a/diffrax/_solver/foster_langevin_srk.py
+++ b/diffrax/_solver/foster_langevin_srk.py
@@ -13,6 +13,7 @@
 
 from .._custom_types import (
     AbstractBrownianIncrement,
+    Args,
     BoolScalarLike,
     DenseInfo,
     RealScalarLike,
@@ -37,7 +38,7 @@
 UnderdampedLangevinArgs = tuple[
     UnderdampedLangevinX,
     UnderdampedLangevinX,
-    Callable[[UnderdampedLangevinX], UnderdampedLangevinX],
+    Callable[[UnderdampedLangevinX, Args], UnderdampedLangevinX],
 ]
 
 
@@ -48,7 +49,7 @@ def _get_args_from_terms(
     PyTree,
     PyTree,
     PyTree,
-    Callable[[UnderdampedLangevinX], UnderdampedLangevinX],
+    Callable[[UnderdampedLangevinX, Args], UnderdampedLangevinX],
 ]:
     drift, diffusion = terms.terms
     if isinstance(drift, WrapTerm):
@@ -255,6 +256,7 @@ def init(
         evaluation of grad_f.
         """
         drift, diffusion = terms.terms
+        del diffusion
         (
             gamma_drift,
             u_drift,
@@ -265,6 +267,7 @@ def init(
 
         h = drift.contr(t0, t1)
         x0, v0 = y0
+        del v0
 
         gamma = broadcast_underdamped_langevin_arg(gamma_drift, x0, "gamma")
         u = broadcast_underdamped_langevin_arg(u_drift, x0, "u")
@@ -287,7 +290,7 @@ def compare_args_fun(arg1, arg2):
         u = jtu.tree_map(compare_args_fun, u, u_diffusion)
 
         try:
-            grad_f_shape = jax.eval_shape(grad_f, x0)
+            grad_f_shape = jax.eval_shape(grad_f, x0, args)
         except ValueError:
             raise RuntimeError(
                 "The function `grad_f` in the Underdamped Langevin term must be"
@@ -300,7 +303,7 @@ def shape_check_fun(_x, _g, _u, _fx):
 
         if not jtu.tree_all(jtu.tree_map(shape_check_fun, x0, gamma, u, grad_f_shape)):
             raise RuntimeError(
-                "The shapes and PyTree structures of x0, gamma, u, and grad_f(x0)"
+                "The shapes and PyTree structures of x0, gamma, u, and grad_f(x0, args)"
                 " must match."
             )
 
@@ -311,7 +314,7 @@ def shape_check_fun(_x, _g, _u, _fx):
 
         coeffs = self._recompute_coeffs(h, gamma, tay_coeffs)
         rho = jtu.tree_map(lambda c, _u: jnp.sqrt(2 * c * _u), gamma, u)
-        prev_f = grad_f(x0) if self._is_fsal else None
+        prev_f = grad_f(x0, args) if self._is_fsal else None
 
         state_out = SolverState(
             gamma=gamma,
@@ -336,6 +339,7 @@ def _compute_step(
         coeffs: _Coeffs,
         rho: UnderdampedLangevinX,
         prev_f: Optional[UnderdampedLangevinX],
+        args: Args,
     ) -> tuple[
         UnderdampedLangevinX,
         UnderdampedLangevinX,
@@ -369,7 +373,6 @@ def step(
     ) -> tuple[
         UnderdampedLangevinTuple, _ErrorEstimate, DenseInfo, SolverState, RESULTS
     ]:
-        del args
         st = solver_state
         drift, diffusion = terms.terms
 
@@ -404,12 +407,12 @@ def step(
             prev_f = st.prev_f
         else:
             prev_f = lax.cond(
-                eqxi.unvmap_any(made_jump), lambda: grad_f(x0), lambda: st.prev_f
+                eqxi.unvmap_any(made_jump), lambda: grad_f(x0, args), lambda: st.prev_f
             )
 
         # The actual step computation, handled by the subclass
         x_out, v_out, f_fsal, error = self._compute_step(
-            h, levy, x0, v0, (gamma, u, grad_f), coeffs, rho, prev_f
+            h, levy, x0, v0, (gamma, u, grad_f), coeffs, rho, prev_f, args
         )
 
         def check_shapes_dtypes(arg, *args):
diff --git a/diffrax/_solver/quicsort.py b/diffrax/_solver/quicsort.py
index 4f21bd6f..dd7c47f6 100644
--- a/diffrax/_solver/quicsort.py
+++ b/diffrax/_solver/quicsort.py
@@ -10,6 +10,7 @@
 
 from .._custom_types import (
     AbstractSpaceTimeTimeLevyArea,
+    Args,
     RealScalarLike,
 )
 from .._local_interpolation import LocalLinearInterpolation
@@ -199,6 +200,7 @@ def _compute_step(
         coeffs: _QUICSORTCoeffs,
         rho: UnderdampedLangevinX,
         prev_f: Optional[UnderdampedLangevinX],
+        args: Args,
     ) -> tuple[UnderdampedLangevinX, UnderdampedLangevinX, None, None]:
         del prev_f
         dtypes = jtu.tree_map(jnp.result_type, x0)
@@ -235,7 +237,7 @@ def _extract_coeffs(coeff, index):
 
         def fn(carry):
             x, _f, _ = carry
-            fx_uh = (f(x) ** ω * uh**ω).ω
+            fx_uh = (f(x, args) ** ω * uh**ω).ω
             return x, _f, fx_uh
 
         def compute_x2(carry):
diff --git a/diffrax/_solver/should.py b/diffrax/_solver/should.py
index caab54d3..4999b9de 100644
--- a/diffrax/_solver/should.py
+++ b/diffrax/_solver/should.py
@@ -6,6 +6,7 @@
 
 from .._custom_types import (
     AbstractSpaceTimeTimeLevyArea,
+    Args,
     RealScalarLike,
 )
 from .._local_interpolation import LocalLinearInterpolation
@@ -198,6 +199,7 @@ def _compute_step(
         coeffs: _ShOULDCoeffs,
         rho: UnderdampedLangevinX,
         prev_f: UnderdampedLangevinX,
+        args: Args,
     ) -> tuple[UnderdampedLangevinX, UnderdampedLangevinX, UnderdampedLangevinX, None]:
         dtypes = jtu.tree_map(jnp.result_type, x0)
         w: UnderdampedLangevinX = jtu.tree_map(jnp.asarray, levy.W, dtypes)
@@ -225,7 +227,7 @@ def _compute_step(
 
         def fn(carry):
             x, _f, _ = carry
-            fx = f(x)
+            fx = f(x, args)
             return x, _f, fx
 
         def compute_x2(carry):
diff --git a/test/test_brownian.py b/test/test_brownian.py
index 3a265019..1acbcf0c 100644
--- a/test/test_brownian.py
+++ b/test/test_brownian.py
@@ -123,7 +123,7 @@ def is_tuple_of_ints(obj):
 def test_statistics(ctr, levy_area, use_levy):
     # Deterministic key for this test; not using getkey()
     key = jr.PRNGKey(5678)
-    num_samples = 60000
+    num_samples = 600000
     keys = jr.split(key, num_samples)
     t0, t1 = 0.0, 5.0
     dt = t1 - t0
@@ -581,7 +581,7 @@ def test_whk_interpolation(tol, spline):
     u = jnp.array(5.7, dtype=jnp.float64)
     bound = 0.0
     rs = jr.uniform(
-        r_key, (100,), dtype=jnp.float64, minval=s + bound, maxval=u - bound
+        r_key, (1000,), dtype=jnp.float64, minval=s + bound, maxval=u - bound
     )
     path = diffrax.VirtualBrownianTree(
         t0=s,
@@ -672,8 +672,8 @@ def eval_paths(t):
     assert jnp.all(_pvals_w > 0.1 / _pvals_w.shape[0])
     assert jnp.all(_pvals_h > 0.1 / _pvals_h.shape[0])
     assert jnp.all(_pvals_k > 0.1 / _pvals_k.shape[0])
-    assert jnp.all(jnp.abs(total_mean_err) < 0.005)
-    assert jnp.all(jnp.abs(total_cov_err) < 0.005)
+    assert jnp.all(jnp.abs(total_mean_err) < 0.01)
+    assert jnp.all(jnp.abs(total_cov_err) < 0.01)
 
 
 def test_levy_area_reverse_time():
diff --git a/test/test_integrate.py b/test/test_integrate.py
index 555d6ade..424146e5 100644
--- a/test/test_integrate.py
+++ b/test/test_integrate.py
@@ -319,7 +319,7 @@ def get_dt_and_controller(level):
         levy_area=None,
         ref_solution=None,
     )
-    assert -0.2 < order - theoretical_order < 0.2
+    assert -0.3 < order - theoretical_order < 0.3
 
 
 # Step size deliberately chosen not to divide the time interval
diff --git a/test/test_sde1.py b/test/test_sde1.py
index b4504872..b50d014f 100644
--- a/test/test_sde1.py
+++ b/test/test_sde1.py
@@ -89,10 +89,9 @@ def get_dt_and_controller(level):
         levy_area=None,
         ref_solution=None,
     )
-    # The upper bound needs to be 0.25, otherwise we fail.
-    # This still preserves a 0.05 buffer between the intervals
-    # corresponding to the different orders.
-    assert -0.2 < order - theoretical_order < 0.25
+    # TODO: this is a pretty wide range to check. Maybe fixable by being better about
+    # the randomness (e.g. average over multiple original seeds)?
+    assert -0.4 < order - theoretical_order < 0.4
 
 
 # Make variables to store the correct solutions in.
diff --git a/test/test_underdamped_langevin.py b/test/test_underdamped_langevin.py
index e945cad5..c0dddec6 100644
--- a/test/test_underdamped_langevin.py
+++ b/test/test_underdamped_langevin.py
@@ -234,7 +234,8 @@ def test_reverse_solve(solver_cls):
 
 # Here we check that if the drift and diffusion term have different arguments,
 # an error is thrown.
-def test_different_args():
+@pytest.mark.parametrize("solver_cls", _only_uld_solvers_cls())
+def test_different_args(solver_cls):
     x0 = (jnp.ones(2), jnp.zeros(2))
     v0 = (jnp.zeros(2), jnp.zeros(2))
     y0 = (x0, v0)
@@ -242,7 +243,7 @@ def test_different_args():
     u1 = (jnp.array([1, 2]), 1)
     g2 = (jnp.array([1, 2]), jnp.array([1, 3]))
     u2 = (jnp.array([1, 2]), jnp.ones((2,)))
-    grad_f = lambda x: x
+    grad_f = lambda x, args: x
 
     w_shape = (
         jax.ShapeDtypeStruct((2,), jnp.float64),
@@ -267,7 +268,7 @@ def test_different_args():
     diffusion_term_b = diffrax.UnderdampedLangevinDiffusionTerm(g1, u2, bm)
     terms_b = diffrax.MultiTerm(drift_term, diffusion_term_b)
 
-    solver = diffrax.ShOULD(0.01)
+    solver = solver_cls(0.01)
     with pytest.raises(Exception):
         diffeqsolve(terms_a, solver, 0, 1, 0.1, y0, args=None)
     diffeqsolve(terms_b, solver, 0, 1, 0.1, y0, args=None)