% The function is used to find the steady state of the model for a chosen gammap. In
% particular, it finds the residual from the specified guess of x(1) that 
% captures how far we are from the steady state for this particular guess
function [resid, rl_rd_spreadp,r_l_subsidyp] = function_r_l(x,betap,xip,deltap,alphap,sigma,sigma_lower_bar,taup,sigmap_zero,sigmap_d,sigmap_c,kappap,pip,f_costs,kappab,Tax1,gammap)

try
rr_rbc= 1/betap-1;
h_rbc = 1;
y_rbc = (alphap/(rr_rbc+deltap))^(alphap/(1-alphap))*h_rbc;
k_rbc = alphap*y_rbc/(rr_rbc+deltap);
    



r_e = (1+kappab)/betap;

l = ((x(1)-(1-deltap))/alphap)^(1/(alphap-1));

%l = k_rbc;



%r_l_subsidyp = 1-(k_rbc^(alphap-1)*alphap+(1-deltap))/x(1);


r_l_subsidyp = 0;

k_r = sigma*l;

k = l*(1-sigma);

invest = deltap*(k+k_r);

h = k/(k+k_r);

h_r = 1 - h;

y = k^(alphap)*h^(1-alphap);

y_r = k_r^alphap*h_r^(1-alphap) + 0.5*k_r - 0.5*(1 + 2*xip)*k_r;

ytot = y + y_r;

e = gammap*l;

d = l-e;

w = (1-alphap)*y/h;



a = 1;


% Old
%Tax1 = 0;
c = (y+y_r)-invest-kappab*e;



lambda_c = (1-betap*kappap)*(c-kappap*c)^(-sigmap_c);

r_d =-(sigmap_zero*d^(-sigmap_d)/lambda_c-1)/betap;


% New
%%d + e + c = w + r_d*d + r_e*e - Tax*ytot;
%%ytot = c + invest + Tax*ytot;

%%d + e + ytot = w + r_d*d + r_e*e + invest;

% r_d = (d + e + ytot - w - r_e*e - invest)/d;
% 
% lambda_c = (sigmap_zero*d^(-sigmap_d))/(1-r_d*betap);
% 
% c = (lambda_c/(1-betap*kappap))^(-1/sigmap_c)/(1-kappap);
% 
% Tax = (ytot - c - invest)/ytot;

% by parametric restriction



rl_rd_spreadp = x(1)- r_d;

    

resid(1) =  -r_e*gammap ... 
            + ( sigma_lower_bar*taup/sqrt(2*pip)...
            *exp(-(  ( f_costs+(x(1)-rl_rd_spreadp)*(1-gammap)-x(1)+xip*sigma_lower_bar )/...
                      (sigma_lower_bar*sqrt(2)*taup)...                
                   )^2 ... 
                 )...
             + 0.5*( x(1)-sigma_lower_bar*xip-(1-gammap)*(x(1)-rl_rd_spreadp)-f_costs )...
                  *(1-erf( ( f_costs+(x(1)-rl_rd_spreadp)*(1-gammap)-x(1)+xip*sigma_lower_bar )/...
                      (sigma_lower_bar*sqrt(2)*taup)... 
                         )...  
                    )...
                );
            

catch
    resid =[1000];
end