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