docs: def eval_components(self, *args, **kwargs):
          """
          Numerically integrate the system of ODEs.
  
          :param args: Ordered arguments for the parameters and independent
            variables
          :param kwargs:  Keyword arguments for the parameters and independent
            variables
          :return:
          """
          bound_arguments = self.__signature__.bind(*args, **kwargs)
          t_like = bound_arguments.arguments[self.independent_vars[0].name]
  
          # System of functions to be integrated
          f = lambda ys, t, *a: [c(t, *(list(ys) + list(a))) for c in self._ncomponents]
          Dfun = lambda ys, t, *a: [[c(t, *(list(ys) + list(a))) for c in row] for row in self._njacobian]
  
          initial_dependent = [self.initial[var] for var in self.dependent_vars]
          t_initial = self.initial[self.independent_vars[0]] # Assuming there's only one
  
          # Check if the time-like data includes the initial value, because integration should start there.
          try:
              t_like[0]
          except (TypeError, IndexError): # Python scalar gives TypeError, numpy scalars IndexError
              t_like = np.array([t_like]) # Allow evaluation at one point.
  
          # The strategy is to split the time axis in a part above and below the
          # initial value, and to integrate those seperately. At the end we rejoin them.
          # np.flip is needed because odeint wants the first point to be t_initial
          # and so t_smaller is a declining series.
          if t_initial in t_like:
              t_bigger = t_like[t_like >= t_initial]
              t_smaller = t_like[t_like <= t_initial][::-1]
          else:
              t_bigger = np.concatenate(
                  (np.array([t_initial]), t_like[t_like > t_initial])
              )
              t_smaller = np.concatenate(
                  (np.array([t_initial]), t_like[t_like < t_initial][::-1])
              )
          # Properly ordered time axis containing t_initial
          t_total = np.concatenate((t_smaller[::-1][:-1], t_bigger))
  
          ans_bigger = odeint(
              f,
              initial_dependent,
              t_bigger,
              args=tuple(
                  bound_arguments.arguments[param.name] for param in self.params),
              Dfun=Dfun,
              *self.lsoda_args, **self.lsoda_kwargs
          )
          ans_smaller = odeint(
              f,
              initial_dependent,
              t_smaller,
              args=tuple(
                  bound_arguments.arguments[param.name] for param in self.params),
              Dfun=Dfun,
              *self.lsoda_args, **self.lsoda_kwargs
          )
  
          ans = np.concatenate((ans_smaller[1:][::-1], ans_bigger))
          if t_initial in t_like:
              # The user also requested to know the value at t_initial, so keep it.
              return ans.T
          else:
              # The user didn't ask for the value at t_initial, so exclude it.
              # (t_total contains all the t-points used for the integration,
              # and so is t_like with t_initial inserted at the right position).
              return ans[t_total != t_initial].T