for t = 1:length(CO2Emissions)

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

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

    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;
end