function [t_out,x_out] = E2depICIpalbo_model_alternation...
    (odefun,tspan,x0,PAR,treatment,cycletime,changemedia,E2in,E2mediaindex,RelTolval,AbsTolval)
tin = tspan(1);
t_out = [];
x_out = [];
while tin < max(tspan)
    tmod = mod(tin,cycletime);
    tcout = 0;
    for i = 1 : numel(treatment)
        if tmod >= tcout && tmod < tcout + treatment{i}.duration * 24
            treatind = i;
        end
        tcout = tcout + treatment{i}.duration * 24;
    end
    
    if ~isempty(t_out)
        x0 = x_out(end,:);
        x0(end) = x0(end);
        x0(E2mediaindex) = treatment{treatind}.E2;
        E2in = treatment{treatind}.E2;
    end
    
    ICI = treatment{treatind}.ICI;
    palbo = treatment{treatind}.palbo;
    abema = treatment{treatind}.abema;
    duration = treatment{treatind}.duration;
    
    tduration = tin:.1:tin + duration * 24;
    ind1 = changemedia * 24 < tduration(end);
    ind2 = changemedia * 24 > tduration(1);
    ind = ind1 & ind2;
    changemedia1 = changemedia(ind);
    [t1,~,x1] = E2depICIpalbo_sim(odefun,PAR,ICI,palbo,abema,tduration,x0,changemedia1,E2in,RelTolval,AbsTolval);
    
    if ~isempty(t_out)
        t_out = [t_out;t1];
    else
        t_out = t1;
    end
    x_out = [x_out;x1];
    tin = tin + treatment{treatind}.duration * 24;
end

end