Derivative of interpolation polynomial is also correct.

scipy_dae/integrate/_dae/radau.py

@@ -699,8 +699,14 @@ def __init__(self, t_old, t, y_old, Q, yp_old, Qp):
     def _call_impl(self, t):
         x = (t - self.t_old) / self.h
         x = np.atleast_1d(x)
-        p = np.tile(x, (self.order + 1, 1))
-        p = np.cumprod(p, axis=0)
+
+        # factors for nterpolation polynomial and its derivative
+        # p = np.tile(x, (self.order + 1, 1))
+        # p = np.cumprod(p, axis=0)
+        c = np.arange(1, self.order + 2)[:, None]
+        p = x**c
+        dp = (c / self.h) * (x**(c - 1))
+
         # Here we don't multiply by h, not a mistake.
         y = np.dot(self.Q, p)
         yp = np.dot(self.Qp, p)
@@ -710,6 +716,14 @@ def _call_impl(self, t):
             y = np.squeeze(y)
             yp = np.squeeze(yp)
 
+        # compute derivative of interpolation polynomial
+        yp = np.dot(self.Q, dp)
+
+        # # compute derivative of interpolation polynomial
+        # yp = np.zeros_like(y)
+        # for i in range(self.order + 1):
+        #     yp += np.outer(self.Q[:, i] * (i + 1), x**i) / self.h
+
         return y, yp