Fsolve 19 state variables with matching dimensions

I am trying to solve (or find the steady-state of) 19 state variables using 19 equations.

The error message I am getting is that for F(3) does not have matching dimensions on left and right side.

Code is provided below.

clear
clc

   %initial conditions
    Fc = [882e3/24,415e3/24,100e3];
    S_I_0 = 40;
    S_S_0 = 125;
    S_ND_0 = 8;
    S_NH_0 = 30;
    S_NI_0 = 2;
    S_NO_0 = 0.5;
    S_O2_0 = 0.1; %0-0.5 g/m3  
    S_ALK_0 = 5;
    
    X_I_0 = 100;
    X_S_0 = 250;
    X_ND_0 = 10;
    X_BH_0 = 1000;
    X_BA_0 = 100;
    X_P_0 = 1;

    X_I_0r = 200;
    X_S_0r = 500;
    X_ND_0r = 20;
    X_BH_0r = 2;
    X_BA_0r = 2;
    X_P_0r = 2;

    init = [S_I_0;S_S_0;S_O2_0;S_NO_0;S_NH_0;S_ND_0;S_ALK_0;
        X_I_0;X_S_0;X_BH_0;X_BA_0;X_P_0;X_ND_0;
        X_I_0r;X_S_0r;X_BH_0r;X_BA_0r;X_P_0r;X_ND_0r
    ];  

    x = fsolve(@(x) ASM1(Fc,x),init)
    SASM1 = x(2) + x(9) + x(5);
    XASM1 = x(10) + x(11);

function F = ASM1(Fc,x)

    S_I = x(1);
    S_S = x(2);
    S_O2 = x(3);
    S_NO = x(4);
    S_NH = x(5);
    S_ND = x(6);
    S_ALK = x(7);

    X_I = x(8);
    X_S = x(9);
    X_BH = x(10);
    X_BA = x(11);
    X_P = x(12);
    X_ND = x(13);

    X_Ir = x(14);
    X_Sr = x(15);
    X_BHr = x(16);
    X_BAr = x(17);
    X_Pr = x(18);
    X_NDr = x(19);
        
        %Flow and Volumes
        V = 5000e3; %m3
        %F_00 = 20*sin(t-1) + 190;
        Fin = 15163e3/24; %m3
        Fr = Fc(1); %m3
        F = Fin + Fr;
        Fw = Fc(2); %m3
        %Fe = F - Fr - Fw;
        Vs = 1017.36e3*3.5*2/3; %m3
        %tao = 24;
        
        %Inlet
        S_S_in = 125;
        S_I_in = 40;
        X_S_in = 250;
        X_I_in = 100;
        S_ND_in = 8;
        X_ND_in = 10;
        S_NH_in = 30;
        S_NO_in = 0.5;
        S_O2_in = 0.1;
%         S_NO_in = 0.5;

        %O2 transfer
        a_0 = 0.0018e-3;
        W = Fc(3);
        C_S = 10; %mg/L
        kLa = a_0*W;
        Qin = kLa*(C_S - S_O2);
        
        %kinetics
        Y_H = 0.67;
        Y_A = 0.24;
        Y_S = 0.63;


        n_g = 0.8; %0.8
        n_h = 0.4;
        f_p = 0.08;
        i_XB = 0.086;
        i_XE = 0.06;

        %% Kinetic Constants
        %K_O2 = 0.2;
        K_S = 20;
        K_NO = 0.5;
        K_OA = 0.4;
        K_OH = 0.2;
        K_NH = 1; %0.05
        K_MAX = 0.34;
        K_X = 0.03;


        mumax_A = 0.8;
        mumax_H = 6;

        b_H = 0.62;
        b_A = 0.132;

        k_h = 3;
        k_a = 0.08;
        %p_7 = X_S * k_7

        mu_O2 = S_O2/(K_OH+S_O2);
        mu_OH = K_OH/(K_OH+S_O2);
        mu_S = S_S/(K_S+S_S);
        mu_NO = S_NO/(K_NO+S_NO);
        mu_NH = S_NH/(K_NH + S_NH);
        mu_SH = (X_S/X_BH)/(K_X+(X_S/X_BH));
        
        
        mu_1 = mumax_H * mu_S * mu_O2 * X_BH;
        mu_2 = mumax_H * mu_S * mu_OH * mu_NO * n_h * X_BH;
        mu_3 = mumax_A * mu_NH * mu_O2 * X_BA;
        mu_4 = b_H * X_BH;
        mu_5 = b_A * X_BA;
        mu_6 = k_a * S_ND * X_BH;
        mu_7 = k_h * mu_SH *(mu_O2 + n_h * mu_OH * mu_NO)*X_BH;
        mu_8 = n_g * X_ND/X_S; %p_7
        
        %% dynamic models
        F(1) = S_I_in*Fin/V - S_I*F/V;
        F(2) = S_S_in*Fin/V - S_S*F/V - mu_1/Y_H - mu_2/Y_H + mu_7;
        F(3) = Qin + S_O2_in*Fin/V - S_O2*F/V - (1-Y_H)/Y_H * mu_1 - (4.57 - Y_A)/Y_A * mu_3;
        F(4) = S_NO_in*Fin/V - S_NO*F/V - (1 - Y_H)/(2.86*Y_H) * mu_2 + mu_3 * Y_A;
        F(5) =  S_NH_in*Fin/V - S_NH*F/V - i_XB* mu_1 - i_XB*mu_2 - (i_XB+1/Y_A) * mu_3 + mu_6;
        F(6) = S_ND_in*Fin/V - S_ND*F/V - mu_6 + mu_8;
        F(7) = - S_ALK*F/V - i_XB/14*mu_1 + ((1-Y_H)/(14*2.86*Y_H)-i_XB/14)*mu_2 - (i_XB/14 + 1/(7*Y_A))*mu_3 + mu_6/14;
        
        F(8) = X_I_in*Fin/V + Fr*X_Ir/V - X_I*F/V;
        F(9) = X_S_in*Fin/V + Fr*X_Sr/V - X_S*F/V + (1-f_p)*(mu_4+mu_5) - mu_7;
        F(10) = Fr*X_BHr/V - X_BH*F/V + mu_1 + mu_2 - mu_4;
        F(11) = Fr*X_BAr/V - X_BA*F/V + mu_3 - mu_5;
        F(12) = Fr*X_Pr/V - X_P*F/V + f_p*(mu_4 + mu_5);
        F(13) = X_ND_in*Fin/V + Fr*X_NDr/V - X_ND*F/V + (i_XB-f_p*i_XE)*(mu_4 + mu_5) - mu_8;
           
        F(14) = F/Vs*X_I - (Fr+Fw)/Vs*X_Ir;
        F(15) = F/Vs*X_S - (Fr+Fw)/Vs*X_Sr;
        F(16) = F/Vs*X_BH - (Fr+Fw)/Vs*X_BHr;
        F(17) = F/Vs*X_BA - (Fr+Fw)/Vs*X_BAr;
        F(18) = F/Vs*X_P - (Fr+Fw)/Vs*X_Pr;
        F(19) = F/Vs*X_ND - (Fr+Fw)/Vs*X_NDr;

end

Answer 1

[Torsten]

Rename F in this assignment:

F = Fin + Fr

F is reserved for the 19-element function vector that you return to fsolve.

[Gabriel McQueen]

Great Answer Torsten! Who wouldda thought it was so simple

[Matt J]

@Gabriel McQueen since Torsten has apparently answered your question, please Accept-click the answer.

[Torsten]

Thank you, Matt.