function [Ku,Kd,Lu,Ld,pu,fval,exitflag]=revise_GLC_steadystate_helper(Ku0,Kd0,Lu0,Ld0,pu0,sigma,varphi, beta, alphau, Au, eta, delta, Ad, alphad, gamma);

    options = optimset('Tolfun',1e-20,'MaxFunEvals',1000000000,'MaxIter',1000000000,'Algorithm','Levenberg-Marquardt');  %设置计算精度：函数值的tolerance  'Algorithm','Levenberg-Marquardt'改变算法
    c=[sigma,varphi, beta, alphau, Au, eta, delta, Ad, alphad, gamma];
    x0=[Ku0,Kd0,Lu0,Ld0,pu0];

    fun=@(x)solveY(x,c);
    [ss,fval,exitflag,output]=fsolve(fun,x0,options);
    Ku=ss(1);
    Kd=ss(2);
    Lu=ss(3);
    Ld=ss(4);
    pu=ss(5);
end



function xss=solveY(x,c)
    Ku=x(1);
    Kd=x(2);
    Lu=x(3);
    Ld=x(4);
    pu=x(5);
    
    sigma=c(1);
    varphi=c(2);
    beta=c(3);
    alphau=c(4);
    Au=c(5);
    eta=c(6); 
    delta=c(7);
    Ad=c(8);
    alphad=c(9);
    gamma=c(10);


    R=1/beta;
    L=Lu+Ld;
    Yu=Au*Ku^alphau*Lu^(1-alphau);
    Iu=delta*Ku;
    Id=delta*Kd;
    w=(1-alphau)*pu*Yu/Lu;
    Yd=w*Ld/(gamma*(1-alphad));
    Yud=(1-gamma)*Yd/pu;
    Bu=(pu*(Yu-Iu)-w*Lu)/(R-1-eta);
    Bd=(Yd-pu*Yud-pu*Id-w*Ld)/(R-1);
    c=Yd;
    T=eta*Bu;
    d=Bu+Bd;

    
    xss(1)=Yd-Ad*(Kd^alphad*Ld^(1-alphad))^gamma*Yud^(1-gamma);
    xss(2)=w-c^sigma*L^varphi;
    xss(3)=R-eta-pu*alphau*Yu/Ku+delta-1;
    xss(4)=R-gamma*alphad*Yd/Kd+delta-1;
    xss(5)=T+c+d-w*L-R*d;
    
    xss(6)=Ku-abs(Ku);
    xss(7)=Kd-abs(Kd);
    xss(8)=Lu-abs(Lu);
    xss(9)=Ld-abs(Ld);
    xss(10)=pu-abs(pu);

end