clc
clear all
tic
%F = 0.001:1:1600;
F = 50:10:1600;
V = 10;
Tk = 128.95 + 273.15;
y0 = [1 4 Tk Tk];
Qk = -4495.7;
for k = 1:length(F)
    rpar = [F(k) Qk];
    estados_est(k,:) = fsolve('vanvusse_EE',y0,[],rpar);
end
F_V = F/V;

figure(2)
plot(F_V,estados_est(:,2))
xlabel('F/V (h^-^1)')
ylabel('Cb_e_e (mol/L)')
toc