In the specific the code is the following:
f = @(m_od, effect) myfunrec(m_od, effect, PARAM, COMPR, RES, i);
[ m_od, effect]=fsolve(f, [PARAM.OP.m, PARAM.OP.eff])
With the function myfunrec:
function F = myfunrec(m_od, effect, PARAM, COMPR, RES, j)
m_ridNom_c=PARAM.OP.m*sqrt(8314/29*293.15)/(RES.rec.des.p_1*10^5);
% Definition of the constant of turbine chocking
k=PARAM.OP.m*sqrt(8314/29*(PARAM.OP.tit+273.15))/(RES.rec.des.p_3*10^5);
m_rC = m_od*sqrt(8314/29*(PARAM.OP.Tin+273.15))/(RES.rec.des.p_1*10^5);
beta_od = m_rC*(PARAM.OP.beta/m_ridNom_c);
p_2od = RES.rec.des.p_1*beta_od;
p_5od = pressDrop(p_2od,PARAM.OP.Dp_rec,1);
p_3od = pressDrop(p_5od,PARAM.OP.Dp_reac,1);
gamma_air = 1.4;
R_air = 8.314/29;
Cp = gamma_air/(gamma_air-1)*R_air;
eta_c = feval(COMPR.eta , m_od , beta_od(end));
T_2od = compressor(beta_od,PARAM.OP.Tin,eta_c);
T_3od = 1/(8314/29)*(k*p_3od*10^5/m_od)^2-273.15;
Q_od = j*RES.rec.des.Q_react;
p_6od = pressDrop(PARAM.OP.Pin,PARAM.OP.Dp_out,0);
p_4od = pressDrop(p_6od,PARAM.OP.Dp_rec,0);
T_4od = expander(p_3od,p_4od,T_3od,0.9,m_od);
T_5od = T_2od + effect*(T_4od-T_2od);
T_6od = T_2od+T_4od-T_5od;
DT_ml = T_4od-T_5od;
Q_rec = m_od*Cp*(T_5od-T_2od);
UA = Q_rec/DT_ml;
F(1) = m_od - Q_od/(Cp*(T_3od-T_2od))
F(2) = UA - RES.rec.des.UA ;
The structures PARAM and RES are contained in the main script