diff --git a/docs/_static/drc.png b/docs/_static/drc.png
index 8915ced2..a86e18d7 100644
Binary files a/docs/_static/drc.png and b/docs/_static/drc.png differ
diff --git a/docs/_static/first-order.png b/docs/_static/first-order.png
index 94b7cff6..9527ce2c 100644
Binary files a/docs/_static/first-order.png and b/docs/_static/first-order.png differ
diff --git a/docs/_static/michaelis-menten-dydt.png b/docs/_static/michaelis-menten-dydt.png
index c0d6e52f..a18d59c8 100644
Binary files a/docs/_static/michaelis-menten-dydt.png and b/docs/_static/michaelis-menten-dydt.png differ
diff --git a/docs/_static/michaelis-menten-tof.png b/docs/_static/michaelis-menten-tof.png
index 2e451ddb..87a4d434 100644
Binary files a/docs/_static/michaelis-menten-tof.png and b/docs/_static/michaelis-menten-tof.png differ
diff --git a/docs/_static/michaelis-menten.png b/docs/_static/michaelis-menten.png
index f7f60028..95e30ff8 100644
Binary files a/docs/_static/michaelis-menten.png and b/docs/_static/michaelis-menten.png differ
diff --git a/docs/_static/simple-first-order.png b/docs/_static/simple-first-order.png
index f77f390d..755b1d45 100644
Binary files a/docs/_static/simple-first-order.png and b/docs/_static/simple-first-order.png differ
diff --git a/overreact/api.py b/overreact/api.py
index 9f544ef1..68af547b 100644
--- a/overreact/api.py
+++ b/overreact/api.py
@@ -503,7 +503,7 @@ def get_drc(
     compounds,
     y0,
     t_span=None,
-    method="BDF",
+    method="Radau",
     qrrho=True,
     scale="l mol-1 s-1",
     temperature=298.15,
diff --git a/overreact/simulate.py b/overreact/simulate.py
index ca99573c..0cb42625 100644
--- a/overreact/simulate.py
+++ b/overreact/simulate.py
@@ -20,7 +20,6 @@
 if _found_jax:
     import jax.numpy as jnp
     from jax import jacfwd
-    from jax import jit
     from jax.config import config
 
     config.update("jax_enable_x64", True)
@@ -28,7 +27,7 @@
     jnp = np
 
 
-def get_y(dydt, y0, t_span=None, method="BDF"):
+def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-3, atol=1e-6):
     """Simulate a reaction scheme from its rate function.
 
     This uses scipy's ``solve_ivp`` under the hood.
@@ -48,6 +47,8 @@ def get_y(dydt, y0, t_span=None, method="BDF"):
         Integration method to use. See `scipy.integrade.solve_ivp` for details.
         Kinetics problems are very often stiff and, as such, "RK45" is
         normally unsuited. "Radau", "BDF" or "LSODA" are good choices.
+    rtol, atol : array-like
+        See `scipy.integrade.solve_ivp` for details.
 
     Returns
     -------
@@ -114,7 +115,16 @@ def get_y(dydt, y0, t_span=None, method="BDF"):
     if hasattr(dydt, "jac"):
         jac = dydt.jac
 
-    res = _solve_ivp(dydt, t_span, y0, method=method, dense_output=True, jac=jac)
+    res = _solve_ivp(
+        dydt,
+        t_span,
+        y0,
+        method=method,
+        dense_output=True,
+        rtol=rtol,
+        atol=atol,
+        jac=jac,
+    )
     y = res.sol
 
     def r(t):
@@ -213,7 +223,7 @@ def _dydt(t, y, k=k_adj, M=M):
         return jnp.dot(A, r)
 
     if _found_jax:
-        _dydt = jit(_dydt)
+        _dydt = _dydt
 
         def _jac(t, y, k=k_adj, M=M):
             # _jac(t, y)[i, j] == d f_i / d y_j