function [Qss, YHss, val] = find_QYH(R, BETA, SIGMA, PHI, CHI, ...
                                     GAMMA, GAMMA_STAR, THETA, EPSILON, ... 
                                     ALPHA, DELTA, XI_P, BF_GDP)
% due to 'fminsearch', input 'x' is unnecessary

% Set the options of the minimization
options = optimset('PlotFcns',@optimplotfval,'PlotFcns',@optimplotx,...
'Display','iter','TolFun',1e-20,'MaxIter',1000000,...
'TolX',1e-20,'MaxFunEvals',1000000);
% Set the initial value of each endogenous variable
x0 = [1 1];
[x1, val]=fminsearch(@(x) F_QYH(x, R, BETA, SIGMA, PHI, CHI, ...
                                GAMMA, GAMMA_STAR, THETA, EPSILON, ...
                                ALPHA, DELTA, XI_P, BF_GDP),x0,options);
                    
% Get the value of each variable and store it
Qss  = x1(1);     %Q  in F_QYH
YHss = x1(2);     %YH in F_QYH

end

