diff --git a/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst b/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst
index 39693a45..17a68eae 100644
--- a/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst
+++ b/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst
@@ -84,4 +84,4 @@ We can see that the reaction went to full completion by checking the final
 concentrations:
 
 >>> y(y.t_max)
-array([0.0000, 1.0000])
+array([0.00, 1.00])
diff --git a/overreact/cli.py b/overreact/cli.py
index c7b23d3e..2fc80e1b 100644
--- a/overreact/cli.py
+++ b/overreact/cli.py
@@ -474,7 +474,9 @@ def _yield_kinetics(self):
                 # TODO(schneiderfelipe): the following is inefficient but probably OK
                 y0[self.model.scheme.compounds.index(name)] = quantity
 
-            y, r = api.get_y(dydt, y0=y0)
+            y, r = api.get_y(
+                dydt, y0=y0, method=self.method, rtol=self.rtol, atol=self.atol
+            )
             conc_table = Table(
                 Column("no", justify="right"),
                 Column("compound", justify="left"),
@@ -492,9 +494,13 @@ def _yield_kinetics(self):
                 )
             yield conc_table
 
-            # TODO(schneiderfelipe): we can get a max time now based on the
-            # changes through time: stop when the graph gets boring.
-            t = np.linspace(y.t_min, y.t_max)
+            eps = 1e-4
+            t_max = y.t_max
+            if np.any(r(y.t_min) > eps):
+                while np.all(r(t_max) < eps):
+                    t_max *= 0.9
+
+            t = np.linspace(y.t_min, t_max)
             if self.plot:
                 for i, name in enumerate(self.model.scheme.compounds):
                     if not core.is_transition_state(name):
diff --git a/overreact/simulate.py b/overreact/simulate.py
index 79eddc4c..ea70fd62 100644
--- a/overreact/simulate.py
+++ b/overreact/simulate.py
@@ -27,7 +27,9 @@
     jnp = np
 
 
-def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11):
+def get_y(
+    dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11, max_time=5 * 60
+):
     """Simulate a reaction scheme from its rate function.
 
     This uses scipy's ``solve_ivp`` under the hood.
@@ -49,6 +51,9 @@ def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11):
         normally unsuited. "Radau", "BDF" or "LSODA" are good choices.
     rtol, atol : array-like
         See `scipy.integrade.solve_ivp` for details.
+    max_time : float, optional
+        If `t_span` is not given, an interval will be estimated, but it can't
+        be larger than this parameter.
 
     Returns
     -------
@@ -71,10 +76,10 @@ def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11):
     used to produce a suitable vector of timepoints for, e.g., plotting:
 
     >>> y.t_min, y.t_max
-    (0.0, 10.0)
+    (0.0, 5.0)
     >>> t = np.linspace(y.t_min, y.t_max)
     >>> t
-    array([ 0.        ,  0.20408163,  ...,  9.79591837, 10.        ])
+    array([..., 0.20408163, ..., 4.89795918, ...])
 
     Both `y` and `r` can be used to check concentrations and rates in any
     point in time. In particular, both are vectorized:
@@ -90,26 +95,27 @@ def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11):
     y0 = np.asanyarray(y0)
 
     if t_span is None:
-        n_halflives = 10.0
+        n_halflives = 5.0  # ensure < 5% remaining material in the worst case
 
         halflife_estimate = 1.0
         if hasattr(dydt, "k"):
             halflife_estimate = (
                 np.max(
                     [
-                        np.max(y0) / 2.0,  # zeroth-order half-life
-                        np.log(2.0),  # first-order half-life
-                        1.0 / np.min(y0[np.nonzero(y0)]),  # second-order half-life
+                        np.max(y0) / 2.0,  # zeroth-order halflife
+                        np.log(2.0),  #  first-order halflife
+                        1.0 / np.min(y0[np.nonzero(y0)]),  # second-order halflife
                     ]
                 )
                 / np.min(dydt.k)
             )
+            logger.info(f"largest halflife guess = {halflife_estimate} s")
 
         t_span = [
             0.0,
-            n_halflives * halflife_estimate,
+            min(n_halflives * halflife_estimate, max_time),
         ]
-        logger.info(f"simulation time span = {t_span} s")
+        logger.info(f"simulation time span   = {t_span} s")
 
     jac = None
     if hasattr(dydt, "jac"):
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 45234358..f82e3db8 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -44,7 +44,7 @@ def test_get_y_propagates_reaction_automatically():
     y, r = simulate.get_y(simulate.get_dydt(scheme, [1.0, 1.0]), y0=y0)
 
     assert y.t_min == 0.0
-    assert y.t_max == 1000.0
+    assert y.t_max == 300.0
     assert y(y.t_min) == pytest.approx(y0)
     assert y(y.t_max) == pytest.approx(
         [1.668212890625, 0.6728515625, 0.341787109375], 9e-5
@@ -86,11 +86,11 @@ def test_get_y_conservation_in_equilibria():
     t = np.linspace(y.t_min, y.t_max, num=100)
 
     assert y.t_min == 0.0
-    assert y.t_max == 10.0
+    assert y.t_max == 5.0
     assert y(y.t_min) == pytest.approx(y0)
-    assert y(y.t_max) == pytest.approx([0.5, 0.5], 3e-5)
+    assert y(y.t_max) == pytest.approx([0.5, 0.5], 5e-5)
     assert r(y.t_min) == pytest.approx([-1, 1])
-    assert r(y.t_max) == pytest.approx([0.0, 0.0], abs=3e-5)
+    assert r(y.t_max) == pytest.approx([0.0, 0.0], abs=5e-5)
 
     assert y.t_min == t[0]
     assert y.t_max == t[-1]