Boundary robustness

src/interpolation_utils.jl

@@ -134,22 +134,30 @@ end
 function du_PCHIP(u, t)
     h = diff(u)
     δ = h ./ diff(t)
+    s = sign.(δ)
 
     function _du(k)
-        if k == 1
-            return ((2 * h[1] + h[2]) * δ[1] - h[1] * δ[2]) / (h[1] + h[2])
-        elseif k == length(t)
-            return ((2 * h[end] + h[end - 1]) * δ[end] - h[end] * δ[end - 1]) /
-                   (h[end] + h[end - 1])
+        sₖ₋₁, sₖ = if k == 1
+            s[1], s[2]
+        elseif k == lastindex(t)
+            s[end - 1], s[end]
+        else
+           s[k - 1], s[k] 
         end
-        sₖ₋₁ = sign(δ[k - 1])
-        sₖ = sign(δ[k])
+
         if sₖ₋₁ == 0 && sₖ == 0
             zero(eltype(δ))
         elseif sₖ₋₁ == sₖ
-            w₁ = 2h[k] + h[k - 1]
-            w₂ = h[k] + 2h[k - 1]
-            δ[k - 1] * δ[k] * (w₁ + w₂) / (w₁ * δ[k] + w₂ * δ[k - 1])
+            if k == 1
+                ((2 * h[1] + h[2]) * δ[1] - h[1] * δ[2]) / (h[1] + h[2])
+            elseif k == lastindex(t)
+                ((2 * h[end] + h[end - 1]) * δ[end] - h[end] * δ[end - 1]) /
+                   (h[end] + h[end - 1])
+            else
+                w₁ = 2h[k] + h[k - 1]
+                w₂ = h[k] + 2h[k - 1]
+                δ[k - 1] * δ[k] * (w₁ + w₂) / (w₁ * δ[k] + w₂ * δ[k - 1])
+            end
         else
             zero(eltype(δ))
         end