diff --git a/climatecritters/model_critters/bistable_melcher.py b/climatecritters/model_critters/bistable_melcher.py
index e95bbe5..35e69c5 100644
--- a/climatecritters/model_critters/bistable_melcher.py
+++ b/climatecritters/model_critters/bistable_melcher.py
@@ -116,84 +116,67 @@ class BistableMelcherModel(CCModel):
 
     Examples
     --------
-    **Generate a synthetic DO record and inspect states**
-
-    .. code-block:: python
-
-        import numpy as np
-        from climatecritters.model_critters.bistable_melcher import (
-            BistableMelcherModel, classify_bistable_states
-        )
-
-        model = BistableMelcherModel(sigma=0.2, gamma=1.5, alpha=-0.4)
-        output = model.integrate(
-            t_span=(0, 599.88), y0=[1.0, 0.0],
-            method='heun_maruyama', dt=0.012,
-            kwargs={'random_seed': 42, 'si': 0.12},
-        )
-
-        # Δb time series as a pyleoclim Series — ready to plot or analyse
-        ts = output.to_pyleo('db')
-        ts.plot()
-
-        # Stadial/interstadial classification is computed automatically
-        states = output.diagnostic_variables['states']   # 1 = cold, 0 = warm
-        print('stadial threshold:',      model.stadial_threshold)
-        print('interstadial threshold:', model.interstadial_threshold)
-
-        # Count DO events (interstadial onsets)
-        n_events = int(np.sum(np.diff(states) == -1))
-        print('number of DO events:', n_events)
-
-    **Inspect all parameters interactively**
-
-    .. code-block:: python
-
-        model.doc('parameters')   # pretty-prints names, current values, descriptions
-
-    **Time-varying forcing** (e.g. a slow CO2 ramp across a glacial cycle)
-
-    .. code-block:: python
-
-        import numpy as np
-        from climatecritters.model_critters.bistable_melcher import BistableMelcherModel
-
-        t_arr = np.linspace(0, 599.88, 5000)
-        gamma_ramp = np.linspace(0.8, 3.2, 5000)   # γ increases over time
-
-        model = BistableMelcherModel(
-            sigma=0.2,
-            gamma=lambda t: float(np.interp(t, t_arr, gamma_ramp)),
-            alpha=0.0,
-        )
-        output = model.integrate(
-            t_span=(t_arr[0], t_arr[-1]), y0=[1.0, 0.0],
-            method='heun_maruyama', dt=0.012,
-            kwargs={'random_seed': 0, 'si': 0.12},
-        )
-        ts = output.to_pyleo('db')
-
-    **Reclassify after adding proxy noise**
-
-    .. code-block:: python
-
-        import numpy as np
-        from climatecritters.model_critters.bistable_melcher import classify_bistable_states
-
-        noisy_db = ts.value + np.random.default_rng(7).normal(0, 0.1, len(ts.value))
-        states_noisy = classify_bistable_states(noisy_db, alpha=-0.4)
-
-    **Sensitivity to calibration constants**
-
-    .. code-block:: python
-
-        model = BistableMelcherModel(b0=0.6, q0=-8.5, q1=11.0, tau=0.95)
-        model.doc('parameters')   # verify the new values are registered
+    Generate a synthetic DO record and inspect states
+
+    ```python
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from climatecritters.model_critters.bistable_melcher import (
+        BistableMelcherModel, classify_bistable_states
+    )
+
+    model = BistableMelcherModel(sigma=0.2, gamma=1.5, alpha=-0.4)
+    output = model.integrate(
+        t_span=(0, 599.88), y0=[1.0, 0.0],
+        method='heun_maruyama', dt=0.012,
+        kwargs={'random_seed': 42, 'si': 0.12},
+    )
+
+    # Δb time series as a pyleoclim Series — ready to plot or analyse
+    ts = output.to_pyleo('db')
+    ts.plot()
+    plt.savefig('docs/reference/figures/BistableMelcher_example.png',
+            dpi=150, bbox_inches='tight')
+
+    # Stadial/interstadial classification is computed automatically
+    states = output.diagnostic_variables['states']   # 1 = cold, 0 = warm
+    print('stadial threshold:',      model.stadial_threshold)
+    print('interstadial threshold:', model.interstadial_threshold)
+
+    # Count DO events (interstadial onsets)
+    n_events = int(np.sum(np.diff(states) == -1))
+    print('number of DO events:', n_events)
+    ```
+
+    Add a time-varying forcing (e.g. a slow CO2 ramp across a glacial cycle):
+
+    ```python
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from climatecritters.model_critters.bistable_melcher import BistableMelcherModel
+
+    t_arr = np.linspace(0, 599.88, 5000)
+    gamma_ramp = np.linspace(0.8, 3.2, 5000)   # γ increases over time
+
+    model = BistableMelcherModel(
+        sigma=0.2,
+        gamma=lambda t: float(np.interp(t, t_arr, gamma_ramp)),
+        alpha=0.0,
+    )
+    output = model.integrate(
+        t_span=(t_arr[0], t_arr[-1]), y0=[1.0, 0.0],
+        method='heun_maruyama', dt=0.012,
+        kwargs={'random_seed': 0, 'si': 0.12},
+    )
+    ts = output.to_pyleo('db')
+    ts.plot()
+    plt.savefig('docs/reference/figures/BistableMelcher_ramp_example.png',
+            dpi=150, bbox_inches='tight')
+    ```
 
     References
     ----------
-    Melcher et al. (2025), Clim. Past, 21, 115-132.
-    https://cp.copernicus.org/articles/21/115/2025/
+    Melcher et al. (2025), Clim. Past, 21, 115-132. https://cp.copernicus.org/articles/21/115/2025/
     """
 
     uses_post_history = True
diff --git a/docs/_quarto.yml b/docs/_quarto.yml
index 9038e81..bd7623f 100644
--- a/docs/_quarto.yml
+++ b/docs/_quarto.yml
@@ -123,6 +123,8 @@ website:
                   href: notebooks/model_demos/g24.ipynb
                 - text: "Stommel"
                   href: notebooks/model_demos/stommel.ipynb
+                - text: "Melcher 2025, Bistable"
+                  href: notebooks/model_demos/bistable_melcher.ipynb
             - section: Attractors
               contents:
                 - text: "Roessler"
@@ -238,6 +240,7 @@ quartodoc:
         - Stocker2003BipolarSeesaw
         - Stocker2003ExtendedSeaIceSeesaw
         - Stommel
+        - BistableMelcherModel
 
     - title: Utils
 
diff --git a/docs/figures/BistableMelcher_example.png b/docs/figures/BistableMelcher_example.png
new file mode 100644
index 0000000..e3f8f17
Binary files /dev/null and b/docs/figures/BistableMelcher_example.png differ
diff --git a/docs/figures/BistableMelcher_ramp_example.png b/docs/figures/BistableMelcher_ramp_example.png
new file mode 100644
index 0000000..9b58d1e
Binary files /dev/null and b/docs/figures/BistableMelcher_ramp_example.png differ
diff --git a/docs/figures/Lorenz96_example.png b/docs/figures/Lorenz96_example.png
index 3c9c0fe..458e28b 100644
Binary files a/docs/figures/Lorenz96_example.png and b/docs/figures/Lorenz96_example.png differ
diff --git a/docs/objects.txt b/docs/objects.txt
index d6cce64..76305c6 100644
--- a/docs/objects.txt
+++ b/docs/objects.txt
@@ -183,6 +183,16 @@ climatecritters.Stocker2003ExtendedSeaIceSeesaw py:class 1 api/Stocker2003Extend
 climatecritters.model_critters.stocker2003_bipolar_seesaw.Stocker2003ExtendedSeaIceSeesaw py:class 1 api/Stocker2003ExtendedSeaIceSeesaw.html#climatecritters.Stocker2003ExtendedSeaIceSeesaw climatecritters.Stocker2003ExtendedSeaIceSeesaw
 climatecritters.Stommel py:class 1 api/Stommel.html#climatecritters.Stommel -
 climatecritters.model_critters.stommel.Stommel py:class 1 api/Stommel.html#climatecritters.Stommel climatecritters.Stommel
+climatecritters.BistableMelcherModel.compute_stability_thresholds py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.compute_stability_thresholds -
+climatecritters.model_critters.bistable_melcher.BistableMelcherModel.compute_stability_thresholds py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.compute_stability_thresholds climatecritters.BistableMelcherModel.compute_stability_thresholds
+climatecritters.BistableMelcherModel.dydt py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.dydt -
+climatecritters.model_critters.bistable_melcher.BistableMelcherModel.dydt py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.dydt climatecritters.BistableMelcherModel.dydt
+climatecritters.BistableMelcherModel.populate_diagnostics_from_history py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.populate_diagnostics_from_history -
+climatecritters.model_critters.bistable_melcher.BistableMelcherModel.populate_diagnostics_from_history py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.populate_diagnostics_from_history climatecritters.BistableMelcherModel.populate_diagnostics_from_history
+climatecritters.BistableMelcherModel.sde_noise py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.sde_noise -
+climatecritters.model_critters.bistable_melcher.BistableMelcherModel.sde_noise py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.sde_noise climatecritters.BistableMelcherModel.sde_noise
+climatecritters.BistableMelcherModel py:class 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel -
+climatecritters.model_critters.bistable_melcher.BistableMelcherModel py:class 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel climatecritters.BistableMelcherModel
 climatecritters.utils.noise.from_series py:function 1 api/utils.noise.from_series.html#climatecritters.utils.noise.from_series -
 climatecritters.utils.noise.from_param py:function 1 api/utils.noise.from_param.html#climatecritters.utils.noise.from_param -
 climatecritters.utils.resample.downsample py:function 1 api/utils.resample.downsample.html#climatecritters.utils.resample.downsample -