for t = 1:length(CO2Emissions)

    NPP = NPP0*(1 + beta*log(B1(t)/B_10));
    B1_dot(t) = alpha_31*B3(t) + alpha_21*B2(t) - NPP + CO2Emissions(t);
    B2_dot = NPP -(alpha_23+alpha_21)*B2(t);
    B2(t+1) = B2(t) + B2_dot;
    B3_dot = alpha_23*B2(t)-alpha_31*B3(t);
    B3(t+1) = B3(t)+ B3_dot;

    for tilde = 1:t
        for i = 1:5
            inre_summa(i) = sum(A(i)*exp(-((t+1)-tilde)/(tau_init(i)*k*sum(B1_dot(1:t)))));
        end
        impulssvar(tilde) = sum(inre_summa)*B1_dot(tilde);
    end

    if t ~= length(CO2Emissions)
        B1(t+1) = B_10 + sum(impulssvar) + sum(A) * (alpha_31*B3(t+1) + alpha_21*B2(t+1) - NPP + CO2Emissions(t+1));
    else
        B1(t+1) = B_10 + sum(impulssvar);
    end

end