code: def chained_parameter_variation(subject, durations, y0, varied_params, default_params=None,
                                  integrate_kwargs=None, x0=None, npoints=1, numpy=None):
      assert len(durations) > 0, 'need at least 1 duration (preferably many)'
      assert npoints > 0, 'need at least 1 point per duration'
      for k, v in varied_params.items():
          if len(v) != len(durations):
              raise ValueError("Mismathced lengths of durations and varied_params")
      if isinstance(subject, ODESys):
          integrate = subject.integrate
          numpy = numpy or subject.numpy
      else:
          integrate = subject
          numpy = numpy or np
      default_params = default_params or {}
      integrate_kwargs = integrate_kwargs or {}
      def _get_idx(cont, idx):
          if isinstance(cont, dict):
              return {k: (v[idx] if hasattr(v, '__len__') and getattr(v, 'ndim', 1) > 0 else v)
                      for k, v in cont.items()}
          else:
              return cont[idx]
      durations = numpy.cumsum(durations)
      for idx_dur in range(len(durations)):
          params = copy.copy(default_params)
          for k, v in varied_params.items():
              params[k] = v[idx_dur]
          if idx_dur == 0:
              if x0 is None:
                  x0 = durations[0]*0
              out = integrate(numpy.linspace(x0, durations[0], npoints + 1), y0, params, **integrate_kwargs)
          else:
              if isinstance(out, Result):
                  out.extend_by_integration(durations[idx_dur], params, npoints=npoints, **integrate_kwargs)
              else:
                  for idx_res, r in enumerate(out):
                      r.extend_by_integration(durations[idx_dur], _get_idx(params, idx_res),
                                              npoints=npoints, **integrate_kwargs)
      return out