How to minimize a variable by varying parameters in a model?

I'm trying to build a model for the dissolution of a molecule in an alkali solution. I started doing it by fixing the parameters I'd like to estimate, buiding the ODE and the iteractions loops that were necessary. I've tried turning the script into a function and making a second script to optimize the STAT, of which is what I want to find the minimum value. The parameters I wish to vary are CS, EA, K0 and pKa. Can anyone guide me into turning the script into something that allows for this minimization?

% Experimental data

T = [67.5, 55, 42.5, 67.5, 67.5, 55, 30, 80, 42.5, 67.5, 55, 42.5, 42.5, 55, 67.5, 42.5, 67.5, 55, 42.5, 42.5, 67.5, 55, 42.5, 67.5, 55, 55, 55, 55];
t = [2.5, 2, 1.5, 1.5, 2.5, 2, 2, 2, 1.5, 2.5, 3, 1.5, 2.5, 2, 2.5, 1.5, 1.5, 1, 2.5, 2.5, 1.5, 2, 2.5, 1.5, 2, 2, 2, 2];
C = [0.2, 0.3, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.2, 0.2, 0.3, 0.2, 0.4, 0.3, 0.4, 0.4, 0.2, 0.3, 0.4, 0.2, 0.4, 0.3, 0.2, 0.2, 0.5, 0.3, 0.3, 0.1]; 
P = [525, 450, 525, 375, 525, 600, 450, 450, 375, 375, 450, 525, 375, 450, 375, 375, 525, 450, 525, 525, 525, 450, 375, 375, 450, 300, 450, 450];  

r = [1034, 1044, 1063, 1052, 1063, 1054, 1044, 1044, 1026, 1026, 1043, 1031, 1051, 1044, 1052, 1051, 1033, 1042, 1063, 1032, 1063, 1044, 1025, 1027, 1070, 1036, 1044, 1014];
xL = [0.085, 0.1118, 0.1706, 0.1222, 0.1673, 0.1444, 0.1118, 0.1099, 0.0596, 0.0609, 0.1056, 0.0852, 0.1219, 0.1125, 0.1217, 0.1228, 0.0834, 0.1077, 0.1632, 0.0857, 0.1689, 0.1113, 0.0604, 0.0603, 0.1783, 0.0750, 0.1122, 0.0395];
pOH = [4.19, 3.40, 3.67, 2.81, 4.09, 4.12, 3.21, 3.59, 2.92, 3.12, 3.45, 3.72, 2.93, 3.45, 2.95, 2.83, 4.00, 3.34, 3.61, 3.76, 3.95, 3.45, 2.91, 3.09, 3.40, 2.44, 3.43, 3.63];
h = [1.57, 2.1, 5.61, 2.4, 5.11, 3.54, 2.07, 1.96, 1.36, 1.35, 1.92, 1.6, 2.38, 2.11, 2.45, 2.45, 1.52, 1.96, 4.92, 1.62, 5.24, 2.05, 1.37, 1.42, 6.28, 1.53, 2.08, 1.24];

% Known parameters

Vag = 0.1; 
rag = 997.04; 
MnNaOH = 40; 
pKw = 14;
R = 8.3144; 
a = 9.55; 
rL = 1365.3;
MML = 960; 
N = 3.82;

% Initial conditions

pOH0 = -log10(C); 
MLI0 = Vag * C .* P; 
Msol0 = Vag * rag + MnNaOH * Vag * C; 

    % Initial specific gravity
    
    nv = length(C); 
    r0 = zeros(1, nv); 

    for i = 1:nv 
    f3_alvo = C(i);
    w = 0.01; 
    erro_C = @(w)(calcular_f3(w, calcular_f2(w, calcular_f1(w))) - f3_alvo)^2; 
    options = optimset('TolX', 1E-6, 'TolFun', 1E-6, 'MaxIter', 500, 'MaxFunEvals', 1000); 
    w_ot = fminsearch(erro_C, w, options); 
    f1_ot = calcular_f1(w_ot); 
    r0(i) = calcular_f2(w_ot, f1_ot); 
    end

V0 = Msol0 ./ r0;     
t0 = zeros(1, nv); 
Ln = zeros(1, nv); 
HnL = zeros(1, nv);

    % Initial guesses

    CS = 0.0001; 
    EA = 50000; 
    K0 = 10000; 
    pKa = 9.5;

    % Vectors

    rPV = zeros(1, nv);
    xLPV = zeros(1, nv);
    MPV = zeros(1, nv);
    VPV = zeros(1, nv);
    NaPV = zeros(1, nv);
    LSV  = zeros(1, nv);
    pOHP = zeros(1, nv); 
    LnP = zeros(1, nv);
    HnLP = zeros(1, nv);

for i = 1:nv  % ODE

    EDO = @(temp, LS) -K0 * exp(-EA/(R*(T(i) + 273.15))) * a * (MLI0(i) - (MML * Msol0(i))/(rL - (MML * (rL/r0(i)) * LS))) * (CS - LS) * ((rL - (MML * (rL/r0(i)) * LS))^2)/(Msol0(i) * (rL - (1 + (MML * (rL/r0(i)))) * LS)); 
    temp = [t0(i) t(i)]; 
    LS = Ln(i) + HnL(i); 
    options = odeset('RelTol',1e-8,'AbsTol',1e-8);
    [temp, LS] = ode45(EDO, temp, LS, options); 
    LSF = LS(end);
    rP = r0(i) + ((rL - r0(i))/r0(i)) * MML * LS(end); 
    xLP = LS(end) * MML / rP; 
    MP = Msol0(i) / (1 - xLP); 
    VP = MP / rP; 
    NaP = C(i) * V0(i) / VP; 
    LSV(i) = LSF;
    rPV(i) = rP;
    xLPV(i) = xLP;
    MPV(i) = MP;
    VPV(i) = VP;
    NaPV(i) = NaP;

end

for i = 1:nv  % Loop pOH <-> LnP

    F = @(z) [z(1)^2 + (N * z(2) - NaPV(i)) * z(1) - 10^pKw; LSV(i) - z(2) * (1 + 10^(pKa + N * (-log10(z(1)) - pKw)))]; 
    CIn = [C(i), CS]; 
    options1 = optimoptions('fsolve', 'Display', 'iter', 'Tolx', 1e-10, 'TolFun', 1e-10);
    Fsol = fsolve(F, CIn, options1);
    pOHP(i) = -log10(Fsol(1)); 
    LnP(i) = Fsol(2);
    HnLP(i) = LSV(i) - Fsol(2); 

end

Er_rP = (r - rPV).^2;
Er_pOH = (pOH - pOHP).^2;
Er_xL = (xL - xLP).^2;
rep_r = [1043.6, 1043.9, 1044.1, 1044.2];
rep_pOH = [3.4, 3.45, 3.45, 3.43];
rep_xL = [11.18, 11.25, 11.13, 11.22];
Stat_r = sum(Er_rP)/var(rep_r);
Stat_pOH = sum(Er_pOH)/var(rep_pOH);
Stat_xL = sum(Er_xL)/var(rep_xL);
STAT = Stat_r + Stat_pOH + Stat_xL;

[Torsten]

As a general guideline, I suggest you study first this example of fitting parameters for monod kinetics:

https://uk.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

Maybe things will get clearer afterwards. At the moment, your code is too weird to make sense of it (at least for me).

[Star Strider]

What are your data? It seems as though ‘t’ is your time vector (and independent variable), and everything else is to be fitted by your ‘EDO’ set of differential equations. First, transpose your data from row vectors to a matrix of column vectors, then write ‘EDO’ to return as many column vectors as there are columns in the data matrix.

I appreciate @Torsten citing to that answer, however a more illuminating one may be Parameter Estimation for a System of Differential Equations.

.