docs: def differentiate(self, coord, edge_order=1, datetime_unit=None):
          """ Differentiate with the second order accurate central
          differences.
  
          .. note::
              This feature is limited to simple cartesian geometry, i.e. coord
              must be one dimensional.
  
          Parameters
          ----------
          coord: str
              The coordinate to be used to compute the gradient.
          edge_order: 1 or 2. Default 1
              N-th order accurate differences at the boundaries.
          datetime_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms',
              'us', 'ns', 'ps', 'fs', 'as'}
              Unit to compute gradient. Only valid for datetime coordinate.
  
          Returns
          -------
          differentiated: Dataset
  
          See also
          --------
          numpy.gradient: corresponding numpy function
          """
          from .variable import Variable
  
          if coord not in self.variables and coord not in self.dims:
              raise ValueError('Coordinate {} does not exist.'.format(coord))
  
          coord_var = self[coord].variable
          if coord_var.ndim != 1:
              raise ValueError('Coordinate {} must be 1 dimensional but is {}'
                               ' dimensional'.format(coord, coord_var.ndim))
  
          dim = coord_var.dims[0]
          if _contains_datetime_like_objects(coord_var):
              if coord_var.dtype.kind in 'mM' and datetime_unit is None:
                  datetime_unit, _ = np.datetime_data(coord_var.dtype)
              elif datetime_unit is None:
                  datetime_unit = 's'  # Default to seconds for cftime objects
              coord_var = coord_var._to_numeric(datetime_unit=datetime_unit)
  
          variables = OrderedDict()
          for k, v in self.variables.items():
              if (k in self.data_vars and dim in v.dims and
                      k not in self.coords):
                  if _contains_datetime_like_objects(v):
                      v = v._to_numeric(datetime_unit=datetime_unit)
                  grad = duck_array_ops.gradient(
                      v.data, coord_var, edge_order=edge_order,
                      axis=v.get_axis_num(dim))
                  variables[k] = Variable(v.dims, grad)
              else:
                  variables[k] = v
          return self._replace_vars_and_dims(variables)