clc
clear all
tic
%F = 0.001:1:1600;
F = 50:10:1600;
V = 10;
y0 = [1 4];
Cain = 5.1;
T  = 135;
rpar(1) = Cain;
rpar(2) = T;

op=optimset('Display','off');

for k = 1:length(F)
    rpar(3) = F(k);
    rpar(1) = 4;
    estados_est(k,:) = fsolve('vanvusse_EE_2eq',y0,op,rpar);

    rpar(1) = 5.1;
    estados_est2(k,:) = fsolve('vanvusse_EE_2eq',y0,op,rpar);
    
    rpar(1) = 6;    
    estados_est3(k,:) = fsolve('vanvusse_EE_2eq',y0,op,rpar);
end
F_V = F/V;

figure(2)
plot(F_V,estados_est(:,2),'-.',F_V,estados_est2(:,2),'-',F_V,estados_est3(:,2),':')
xlabel('F/V (h^-^1)')
ylabel('Cb_e_e (mol/L)')
legend('Ca_{in}=4.0 mol/L','Ca_{in}=5.1 mol/L','Ca_{in}=6.0 mol/L');

rpar(1) = Cain;
for k = 1:length(F)
    rpar(3) = F(k);
    rpar(2) = 115;
    estados_est(k,:) = fsolve('vanvusse_EE_2eq',y0,op,rpar);

    rpar(2) = 135;
    estados_est2(k,:) = fsolve('vanvusse_EE_2eq',y0,op,rpar);
    
    rpar(2) = 155;    
    estados_est3(k,:) = fsolve('vanvusse_EE_2eq',y0,op,rpar);
end
F_V = F/V;

figure(3)
plot(F_V,estados_est(:,2),'-.',F_V,estados_est2(:,2),'-',F_V,estados_est3(:,2),':')
xlabel('F/V (h^-^1)')
ylabel('Cb_e_e (mol/L)')
legend('T=115 ^oC','T=135 ^oC','T=155 ^oC');
toc
