solve_interp.m

function [x, r1] = solve_interp(m0, CHCl, m, H, KW, KF, KS, ST, FT, f_guess, at_guess)
    global m0 CHCl m H KW KF KS ST FT
    options = optimset('TolX',1e-32,'TolFun',1e-32,'MaxIter', 1000, 'MaxFunEvals', 5000,)
    x0 = [f_guess, at_guess];
    
    % Use anonymous function to pass constants to interp1
    fun = @(x) interp1(x, m0, CHCl, m, H, KW, KF, KS, ST, FT);
    
    [x, fval] = fsolve(@interp1, x0, options);
    r1 = fval;
end

function r1 = interp1(x)%, m0, CHCl, m, H, KW, KF, KS, ST, FT)
    global m0 CHCl m H KW KF KS ST FT
    f    = x(1); 
    AT   = x(2);
    Z    = 1 + ST/KS; 
    HSO4 = m0.*ST./(1+(Z*KS)./(f*H)); 
    HF   = m0.*FT./(1+KF./(f*H));
    r1   = m0*AT + HSO4 + HF - m.*CHCl + (m0 + m).*((f*H./Z) - (Z.*KW./(f*H)));
end