Problem with initial guess of a nonlinear system using fsolve

Dear all, actually I’m dealing with a problem involving the solution of a nonlinear system. I’m recicling the lines I wrote above for the solution of an implicit equation so:

f=@(x) myfun(x,k1,k2...) With k constants. [x]=fsolve(f,x0)

With myfun: Function F=myfun(x,k1,k2,...)

..some passages..

F(1)=implicit eq.

It works fine. If I use the same form to solve a system:

f=@(x,y) myfun(x,y,k1,k2...) With k constants.

x0y0=[x0,y0] [x,y]=fsolve(f,x0y0)

With myfun:

Function F=myfun(x,y,k1,k2,...)

..some passages..

F(1)=funct(x,y) F(2)=funct(x,y)

I get an error, and it is due to the fact that the vector of inital guess gives problem. I also tried writing fsolve(f,x0,y0) but only x0 is read as initial guess, while y0 is read as option (a structure). Could you suggest me how to fix this problem? Thanks in advance, Niccolò

[niccolò guazzi]

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