global  r0 r1 r2 v1 v2 k1 k2   

colors = {'#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd','#8c564b','#e377c2'};

opt = optimoptions('fsolve','FunctionTolerance', 1e-32,'OptimalityTolerance',1e-32, 'MaxFunctionEvaluation', 10e4, 'MaxIterations', 10e3 ,'Display','none');

s0_range = linspace(10,2500, 100);

PC1 =[];F0_s=[];

    r0 = 0e12;    r1 =  7.5e12; r2 =  7.5e12;
    v1 = 12.5e-9; v2 = 12.5e-9; k1 = 75;   k2 = 75;  

    for i = 1:length(s0_range)
        
        s0 = s0_range(i);
        x0 = [0,0,s0];
        
        [x, fval] = fsolve(@Variable,x0,opt, s0);
        s1=x(1); s2=x(2);sF=x(3);
        
        s0_ =  1e-3 * ((s0+sF)/2);
        s1_ =  1e-3 * ((sF+s1)/2);
        s2_ =  1e-3 * ((sF+s2)/2);
        
        r0_s = r0 *((0.685*s0_^4) - (1.01411*s0_^3) + (0.9513*s0_^2) + (0.1364*s0_) + 0.3396);
        r1_s = r1 *((0.685*s1_^4) - (1.01411*s1_^3) + (0.9513*s1_^2) + (0.1364*s1_) + 0.3396);
        r2_s = r2 *((0.685*s2_^4) - (1.01411*s2_^3) + (0.9513*s2_^2) + (0.1364*s2_) + 0.3396);
        
        F1 = (v1*s1)/(k1+s1) ;
        F2 = (v2*s2)/(k2+s2);
        F0 = F1+F2;
        
        PC1 = [PC1 F1/F0]; F0_s = [F0_s F0];
    end
  
plot(F0_s*1e9, PC1*100, 'LineWidth', 1.5);
xlabel('Loading rate (nmol/s)'); ylabel('PC1 (%)');

function F = Variable(x, s0)
s1=x(1); s2=x(2); sF=x(3);

global  r0 r1 r2 v1 v2 k1 k2

sF=s0;

s0_ = 1e-3 * ((s0+sF)/2);
s1_ = 1e-3 * ((sF+s1)/2);
s2_ = 1e-3 * ((sF+s2)/2);

r0_s = r0 *((0.685*s0_^4) - (1.01411*s0_^3) + (0.9513*s0_^2) + (0.1364*s0_) + 0.3396);
r1_s = r1 *((0.685*s1_^4) - (1.01411*s1_^3) + (0.9513*s1_^2) + (0.1364*s1_) + 0.3396);
r2_s = r2 *((0.685*s2_^4) - (1.01411*s2_^3) + (0.9513*s2_^2) + (0.1364*s2_) + 0.3396);

sF=s0;

U1_flow = s0*( (r2_s*(s0-s1) + r0_s*(s2-s1))) / (r0_s*(r1_s*r2_s) + r1_s*r2_s );
U1_sink = (v1*s1)/(k1+s1) ;

U2_flow = s0*( (r1_s*(s0-s2) + r0_s*(s1-s2))) / (r0_s*(r1_s*r2_s) + r1_s*r2_s ) ;
U2_sink = (v2*s2)/(k2+s2);


F(1) = 1e9 * (U1_flow - U1_sink);
F(2) = 1e9 * (U2_flow - U2_sink);
F(3) = 1e9 * (s0-sF) ;

end