diff --git a/compare_dry_density.py b/compare_dry_density.py
new file mode 100644
index 00000000..a0ecd617
--- /dev/null
+++ b/compare_dry_density.py
@@ -0,0 +1,35 @@
+# %%
+import numpy as np
+import matplotlib.pyplot as plt
+
+# %%
+from PyMPDATA_examples.Shipway_and_Hill_2012 import si, settings
+
+# %%
+RHOD_VERTVELO = 3*si.m/si.s * si.kg/si.m**3
+T_MAX = 15*si.minutes
+P0 = 1000*si.hPa
+Z_MAX = 3200*si.m
+nr = 32
+dt = 1.25*si.s
+dz = 25*si.m
+R_MIN = 1 * si.um
+R_MAX = 20.2 * si.um
+apprx_drhod_dz = False
+
+settings_true = settings.Settings(
+                    rhod_w_const=RHOD_VERTVELO, nr=nr, dt=dt, dz=dz, t_max=T_MAX, 
+                    r_min=R_MIN, r_max=R_MAX, p0=P0, z_max=Z_MAX, apprx_drhod_dz=apprx_drhod_dz
+                )
+settings_false = settings.Settings(
+                    rhod_w_const=RHOD_VERTVELO, nr=nr, dt=dt, dz=dz, t_max=T_MAX, 
+                    r_min=R_MIN, r_max=R_MAX, p0=P0, z_max=Z_MAX, apprx_drhod_dz=apprx_drhod_dz, use_max_step=False,
+                )
+# %%
+z = z_points = np.linspace(0.0, 2 * settings_true.nz + 1)
+plt.plot(settings_false.rhod(z) - settings_true.rhod(z), z)
+plt.legend()
+plt.ylabel("z / m")
+plt.xlabel("difference in dry density compared to using max_step / kg/m$^{-3}$")
+plt.show()
+# %%
diff --git a/examples/PyMPDATA_examples/Shipway_and_Hill_2012/settings.py b/examples/PyMPDATA_examples/Shipway_and_Hill_2012/settings.py
index 1490cfcb..e5d0c3d1 100644
--- a/examples/PyMPDATA_examples/Shipway_and_Hill_2012/settings.py
+++ b/examples/PyMPDATA_examples/Shipway_and_Hill_2012/settings.py
@@ -27,6 +27,7 @@ def __init__(
         ksi_1: float = default_ksi_1.to_base_units().magnitude,
         z_max: float = 3000 * si.metres,
         apprx_drhod_dz: bool = True,
+        use_max_step: bool = True,
     ):
         self.dt = dt
         self.dz = dz
@@ -58,11 +59,16 @@ def drhod_dz(z, rhod):
             return drhod_dz
 
         z_points = np.arange(0, self.z_max + self.dz / 2, self.dz / 2)
+        if use_max_step:
+            max_step=self.dz/2
+        else:
+            max_step = np.inf  # (same as default in solve_ivp docs)
         rhod_solution = solve_ivp(
             fun=drhod_dz,
             t_span=(0, self.z_max),
             y0=np.asarray((self.rhod0,)),
             t_eval=z_points,
+            max_step=max_step,
         )
         assert rhod_solution.success