diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 5f686539b..466e95b08 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -75,24 +75,9 @@ class Function(pybamm.Symbol):
         if variable == self:
             return pybamm.Scalar(1)
         else:
-            children = self.orphans
-            partial_derivatives: list[None | pybamm.Symbol] = [None] * len(children)
-            for i, child in enumerate(self.children):
-                # if variable appears in the function, differentiate
-                # function, and apply chain rule
-                if variable in child.pre_order():
-                    partial_derivatives[i] = self._function_diff(
-                        children, i
-                    ) * child.diff(variable)
-
-            # remove None entries
-            partial_derivatives = [x for x in partial_derivatives if x is not None]
-
-            derivative = sum(partial_derivatives)
-            if derivative == 0:
-                return pybamm.Scalar(0)
-
-            return derivative
+            raise NotImplementedError(
+                "Derivative of base Function class is not implemented"
+            )
 
     def _function_diff(self, children: Sequence[pybamm.Symbol], idx: float):
         """
diff --git a/pybamm/expression_tree/interpolant.py b/pybamm/expression_tree/interpolant.py
index 84c982342..97074a817 100644
--- a/pybamm/expression_tree/interpolant.py
+++ b/pybamm/expression_tree/interpolant.py
@@ -197,6 +197,47 @@ class Interpolant(pybamm.Function):
         self.interpolator = interpolator
         self.extrapolate = extrapolate
 
+    def diff(self, variable):
+        """Compute the derivative of the interpolant with respect to a variable.
+        
+        For 1D interpolants with a single child, the derivative can be computed
+        for cubic interpolators (using scipy's derivative method). For linear and
+        pchip interpolators, as well as multi-child interpolants, differentiation
+        is not implemented.
+        """
+        # If the variable is not one of our children, return 0
+        if variable not in self.orphans:
+            return pybamm.Scalar(0)
+        
+        # Check if we have more than one child - differentiation not implemented
+        if len(self.children) > 1:
+            raise NotImplementedError(
+                "differentiation not implemented for functions with more than one child"
+            )
+        
+        # For 1D interpolants with single child, try to compute derivative
+        if self.dimension == 1:
+            try:
+                # Try to use the scipy function's derivative method
+                # CubicSpline has derivative(), interp1d and PchipInterpolator don't
+                deriv_func = self.function.derivative()
+                return pybamm.Function(
+                    deriv_func,
+                    *self.children,
+                    derivative=self.derivative,
+                    differentiated_function=self.function,
+                )
+            except (AttributeError, NotImplementedError):
+                # If derivative method not available, raise NotImplementedError
+                raise NotImplementedError(
+                    "Differentiation not implemented for this interpolator type"
+                )
+        
+        # For 2D and 3D interpolants, not implemented
+        raise NotImplementedError(
+            "Differentiation not implemented for %dD interpolation" % self.dimension
+        )
+
     @classmethod
     def _from_json(cls, snippet: dict):
         """Create an Interpolant object from JSON data"""
@@ -323,6 +364,7 @@ class Interpolant(pybamm.Function):
             "y": self.y.tolist(),
             "interpolator": self.interpolator,
             "extrapolate": self.extrapolate,
+            "_num_derivatives": 0,
         }
 
         return json_dict