In Fritsch and Carlson's paper on monotone interpolation, they identify numerous conditions under which a cubic Hermite interpolator will be monotone. For example: On the subinterval $[t_i, t_{i+1}]$ where $p(t_i) = y_i$ and $p(t_{i+1}) = y_{i+1}$, if we choose endpoint derivatives $d_i = d_{i+1} = 0$, then the interpolator is monotone. There are other conditions. Let $h_i := t_{i+1}-t_i$ and $\Delta_i := (y_{i+1} - y_{i})/h_i$. The condition $$ d_i + d_{i+1} - 2\Delta_i \le 0 $$ is sufficient to ensure monotonicity. If $y_i' + y_{i+1}' - 2\Delta_i > 0$ then we could either have $x^{*} := \frac{1}{3} \frac{2d_{i} + d_{i+1} - 3\Delta_i}{d_{i} + d_{i+1} - 2 \Delta_i}$ not in the interval $[0, 1]$, or if it is in the interval $[0,1]$, then $\mathrm{sgn}(p(x^{*})) = \mathrm{sgn}(\Delta_i)$ will again guarantee monotonicity.
This is somewhat complicated, but there is quite a bit of freedom on a single interval. However, $d_i$ is constrained by data to both its right and left. It's clear that if $\Delta_{i-1}$ and $\Delta_i$ do not have the same sign, that $d_i = 0$ is the only solution. However, when $\Delta_{i-1}$ and $\Delta_i$ have the same sign, we need to find a choice of $d_i$ that satisfies the constraints listed above on both the intervals $[t_{i-1}, t_i]$ and $[t_i, t_{i+1}]$. Matlab solves this problem by setting $d_i$ to the "weighted harmonic mean" \begin{align} \frac{w_1 + w_2}{d_i} = \frac{w_1}{\Delta_{i-1}} + \frac{w_2}{\Delta_i}, w_1 :=2h_k + h_{k-1}, w_2 := h_k + 2h_{k-1} \end{align} I am unfamiliar with the properties of weighted harmonic means, and hence I cannot see why this choice for the derivatives satisfied the constraints for monotonicity listed in Fritsch and Carlson's paper.
How do we prove that the weighted harmonic mean leads to a monotone interpolator?
Partial solution: Using @MPIchael's suggestion, I was at least able to split of some cases. For example, if $d_{i+1}$ must be zero, then $$ d_i - 2\Delta_i \le \max(\Delta_i, \Delta_{i-1}) - 2\Delta_i. $$ And if $\max(\Delta_i, \Delta_{i-1}) = \Delta_i$, we satisfy $d_i + d_{i+1} - 2\Delta_i \le 0$. Otherwise, use $d_i \le \sqrt{\Delta_i \Delta_{i-1}}$, whereupon $$ d_i - 2\Delta_i \le \Delta_i(\sqrt{\Delta_{i-1}/\Delta_i }-2) $$ which is non-positive whenever $\Delta_{i-1} \le 4\Delta_i$. I suspect to get the more general case, I need the inequality $H_{w}(\Delta_i, \Delta_{i-1}) \le 2\min(\Delta_i, \Delta_{i-1})$, which I know to be true for the standard harmonic mean, but I can't find a proof for the weighted harmonic mean.
