
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% STEADY STATE FOR RAMSEY LAGRANGE MULTIPLIERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% In this file I:
% 1. set the SS eqs for Ramsey FOCs
% 2. set a level for instruments and then obtain the private variables block Steady States (find_ss(tools))
% 3. Substitute the results in the SS Eqs for the Ramsey FOCs

% 4. Set the eqs. in matrix form and solve for the associated SS of the LM
% 5. calculate the error of approximation for the instruments level set initially

function [LMOut, errorNorm] = LMultStaticCoop(tools)

    % Initial setup
    nVars = 54; %Total number of variables (included instruments)
    utilVars = 6; %Number of variables that enter the utility

    %Variables (SS)
    syms Cc Ca Cb Hc Ha Hb R Rdc Rda Rdb Bc Ba Bb Yc Ya Yb Kc Ka Kb Ic Ia Ib...
        Qc Qa Qb Rkc Rka Rkb Nc Na Nb Fa Fb Dc Rba Rbb Lambdac Lambdaa Lambdab...
        jc ja jb muc mua mub tauc taua taub
    syms Ac Aa Ab xic xia xib epsAc epsAa epsAb epsxic epsxia epsxib
    % Steady States
    %syms Dc_bar Fa_bar Fb_bar 

    %LM
        syms Lw1 Lw2 Lw3 Lw4 Lw5 Lw6 Lw7 Lw8 Lw9 Lw10 Lw11 Lw12 Lw13 Lw14 Lw15 Lw16 Lw17 Lw18 Lw19 Lw20 Lw21 Lw22 Lw23 Lw24...
        Lw25 Lw26 Lw27 Lw28 Lw29 Lw30 Lw31 Lw32 Lw33 Lw34 Lw35 Lw36 Lw37 Lw38 Lw39 Lw40 Lw41 Lw42 Lw43 Lw44 Lw45 Lw46 Lw47...
        Lw48 Lw49 Lw50 Lw51 

    % Parameters
    syms nc na nb betta thetta sigmma pssi alphha delta zetta eta etaD deltaT psiDa psiDb kappac kappaa kappab kappafa kappafb...
        rhoA rhok stdepsA stdxik phi_tau phi_cr phih 


    %Static FOCs of Ramsey Planner [Cooperative] 
    
    focW1_ss = Lw42*nc + nc/Cc^sigmma + (Lw49*betta*sigmma)/(Cc*(Cc/Cc)^sigmma) + (Lw7*sigmma*exp(Ac)*(Kc*exp(xic))^alphha*(alphha - 1))/(Cc^(sigmma + 1)*Hc^alphha); 
    focW2_ss = Lw42*na - Lw10 + na/Ca^sigmma + (Lw50*betta*sigmma)/(Ca*(Ca/Ca)^sigmma) + (Lw8*sigmma*exp(Aa)*(Ka*exp(xia))^alphha*(alphha - 1))/(Ca^(sigmma + 1)*Ha^alphha); 
    focW3_ss = Lw42*nb - Lw11 + nb/Cb^sigmma + (Lw51*betta*sigmma)/(Cb*(Cb/Cb)^sigmma) + (Lw9*sigmma*exp(Ab)*(Kb*exp(xib))^alphha*(alphha - 1))/(Cb^(sigmma + 1)*Hb^alphha); 
    focW4_ss = (Kc^(alphha - 1)*Lw18*alphha*exp(Ac)*exp(xic)^alphha*(alphha - 1)*(tauc - 1))/(Hc^alphha*Qc) - Hc^pssi*nc - (Kc^alphha*Lw38*exp(Ac)*exp(xic)^alphha*(alphha - 1))/Hc^alphha - Lw7*(Hc^(pssi - 1)*pssi - (alphha*exp(Ac)*(Kc*exp(xic))^alphha*(alphha - 1))/(Cc^sigmma*Hc^(alphha + 1))); 
    focW5_ss = (Ka^(alphha - 1)*Lw19*alphha*exp(Aa)*exp(xia)^alphha*(alphha - 1)*(taua - 1))/(Ha^alphha*Qa) - Ha^pssi*na - (Lw10*exp(Aa)*(alphha - 1)*((Ka*exp(xia))^alphha - alphha*(Ka*exp(xia))^alphha + Ka^alphha*alphha*taua*exp(xia)^alphha))/Ha^alphha - (Ka^alphha*Lw39*exp(Aa)*exp(xia)^alphha*(alphha - 1))/Ha^alphha - Lw8*(Ha^(pssi - 1)*pssi - (alphha*exp(Aa)*(Ka*exp(xia))^alphha*(alphha - 1))/(Ca^sigmma*Ha^(alphha + 1))); 
    focW6_ss = (Kb^(alphha - 1)*Lw20*alphha*exp(Ab)*exp(xib)^alphha*(alphha - 1)*(taub - 1))/(Hb^alphha*Qb) - Hb^pssi*nb - (Lw11*exp(Ab)*(alphha - 1)*((Kb*exp(xib))^alphha - alphha*(Kb*exp(xib))^alphha + Kb^alphha*alphha*taub*exp(xib)^alphha))/Hb^alphha - (Kb^alphha*Lw40*exp(Ab)*exp(xib)^alphha*(alphha - 1))/Hb^alphha - Lw9*(Hb^(pssi - 1)*pssi - (alphha*exp(Ab)*(Kb*exp(xib))^alphha*(alphha - 1))/(Cb^sigmma*Hb^(alphha + 1))); 
    focW7_ss = betta*(Ba*Lw10 + Bb*Lw11) - Lambdac*Lw1 - Lambdab*Lw3 - Lambdaa*Lw2; 
    focW8_ss = Lambdac*Lw26*(jc*thetta - thetta + 1) - Lambdac*Lw21*(jc*thetta - thetta + 1) - Lambdac*Lw24*(jc*thetta - thetta + 1) - Lambdac*Lw25*(jc*thetta - thetta + 1) - Lambdac*Lw4 - Dc*Lw29*betta*thetta; 
    focW9_ss = (Lambdaa*Lw27*(psiDa - 1)*(ja*thetta - thetta + 1))/psiDa - (Lambdaa*Lw22*(psiDa - 1)*(ja*thetta - thetta + 1))/psiDa - Lambdaa*Lw5 + Fa*betta*thetta*(Lw10 - Lw30)*(psiDa - 1); 
    focW10_ss = (Lambdab*Lw28*(psiDb - 1)*(jb*thetta - thetta + 1))/psiDb - (Lambdab*Lw23*(psiDb - 1)*(jb*thetta - thetta + 1))/psiDb - Lambdab*Lw6 + Fb*betta*thetta*(Lw11 - Lw31)*(psiDb - 1); 
    focW11_ss = Lw1*eta + Lw41*nc + Bc*Lw42*eta*nc; 
    focW12_ss = Lw2*eta - Lw10*(Ba*eta + 1) + Lw41*na + Lw10*R*betta + Ba*Lw42*eta*na; 
    focW13_ss = Lw3*eta - Lw11*(Bb*eta + 1) + Lw41*nb + Lw11*R*betta + Bb*Lw42*eta*nb; 
    focW14_ss = - Lw38 - Lw42*nc; 
    focW15_ss = - Lw39 - Lw42*na; 
    focW16_ss = - Lw40 - Lw42*nb; 
    focW17_ss = Lw43*betta*rhoA - Lw43 + Hc^(1 - alphha)*Kc^alphha*Lw38*exp(Ac)*exp(xic)^alphha - (Lw7*exp(Ac)*(Kc*exp(xic))^alphha*(alphha - 1))/(Cc^sigmma*Hc^alphha) - (Hc^(1 - alphha)*Kc^(alphha - 1)*Lw18*alphha*exp(Ac)*exp(xic)^alphha*(tauc - 1))/Qc; 
    focW18_ss = Lw44*betta*rhoA - Lw10*((Ha*exp(Aa)*(Ka*exp(xia))^alphha*(alphha - 1))/Ha^alphha - Ha^(1 - alphha)*Ka*Ka^(alphha - 1)*alphha*taua*exp(Aa)*exp(xia)^alphha) - Lw44 + Ha^(1 - alphha)*Ka^alphha*Lw39*exp(Aa)*exp(xia)^alphha - (Lw8*exp(Aa)*(Ka*exp(xia))^alphha*(alphha - 1))/(Ca^sigmma*Ha^alphha) - (Ha^(1 - alphha)*Ka^(alphha - 1)*Lw19*alphha*exp(Aa)*exp(xia)^alphha*(taua - 1))/Qa; 
    focW19_ss = Lw45*betta*rhoA - Lw11*((Hb*exp(Ab)*(Kb*exp(xib))^alphha*(alphha - 1))/Hb^alphha - Hb^(1 - alphha)*Kb*Kb^(alphha - 1)*alphha*taub*exp(Ab)*exp(xib)^alphha) - Lw45 + Hb^(1 - alphha)*Kb^alphha*Lw40*exp(Ab)*exp(xib)^alphha - (Lw9*exp(Ab)*(Kb*exp(xib))^alphha*(alphha - 1))/(Cb^sigmma*Hb^alphha) - (Hb^(1 - alphha)*Kb^(alphha - 1)*Lw20*alphha*exp(Ab)*exp(xib)^alphha*(taub - 1))/Qb; 
    focW20_ss = Lw46*betta*rhok - Lw46 - (Lw18*(Qc*exp(xic)*(delta - 1) + Hc^(1 - alphha)*Kc^(alphha - 1)*alphha^2*exp(Ac)*exp(xic)*exp(xic)^(alphha - 1)*(tauc - 1)))/Qc - Kc*Lw15*exp(xic)*(delta - 1) + Hc^(1 - alphha)*Kc^alphha*Lw38*alphha*exp(Ac)*exp(xic)*exp(xic)^(alphha - 1) - (Kc*Lw7*alphha*exp(Ac)*exp(xic)*(Kc*exp(xic))^(alphha - 1)*(alphha - 1))/(Cc^sigmma*Hc^alphha); 
    focW21_ss = Lw10*(Ha^(1 - alphha)*Ka*Ka^(alphha - 1)*alphha^2*taua*exp(Aa + xia)*exp(xia)^(alphha - 1) - (Ha*Ka*alphha*exp(Aa + xia)*(Ka*exp(xia))^(alphha - 1)*(alphha - 1))/Ha^alphha) - Lw47 + Lw47*betta*rhok - (Lw19*(Qa*exp(xia)*(delta - 1) + Ha^(1 - alphha)*Ka^(alphha - 1)*alphha^2*exp(Aa + xia)*exp(xia)^(alphha - 1)*(taua - 1)))/Qa - Ka*Lw16*exp(xia)*(delta - 1) + Ha^(1 - alphha)*Ka^alphha*Lw39*alphha*exp(Aa + xia)*exp(xia)^(alphha - 1) - (Ka*Lw8*alphha*exp(Aa + xia)*(Ka*exp(xia))^(alphha - 1)*(alphha - 1))/(Ca^sigmma*Ha^alphha); 
    focW22_ss = Lw11*(Hb^(1 - alphha)*Kb*Kb^(alphha - 1)*alphha^2*taub*exp(Ab + xib)*exp(xib)^(alphha - 1) - (Hb*Kb*alphha*exp(Ab + xib)*(Kb*exp(xib))^(alphha - 1)*(alphha - 1))/Hb^alphha) - Lw48 + Lw48*betta*rhok - (Lw20*(Qb*exp(xib)*(delta - 1) + Hb^(1 - alphha)*Kb^(alphha - 1)*alphha^2*exp(Ab + xib)*exp(xib)^(alphha - 1)*(taub - 1)))/Qb - Kb*Lw17*exp(xib)*(delta - 1) + Hb^(1 - alphha)*Kb^alphha*Lw40*alphha*exp(Ab + xib)*exp(xib)^(alphha - 1) - (Kb*Lw9*alphha*exp(Ab + xib)*(Kb*exp(xib))^(alphha - 1)*(alphha - 1))/(Cb^sigmma*Hb^alphha); 
    focW23_ss = Lw32*Qc*kappac - betta*(Lw15*exp(xic)*(delta - 1) - Lw29*(Qc*deltaT + Qc*Rkc*thetta) - Hc^(1 - alphha)*Kc^(alphha - 1)*Lw38*alphha*exp(Ac)*exp(xic)^alphha + (Lw7*alphha*exp(Ac)*exp(xic)*(Kc*exp(xic))^(alphha - 1)*(alphha - 1))/(Cc^sigmma*Hc^alphha) + (Hc^(1 - alphha)*Kc^(alphha - 2)*Lw18*alphha*exp(Ac)*exp(xic)^alphha*(alphha - 1)*(tauc - 1))/Qc) - Lw35*Qc - Lw15; 
    focW24_ss = Lw33*Qa*kappaa - betta*(Lw10*(Qa*deltaT + Qa*Rka*(thetta - 1) - Ha^(1 - alphha)*Ka^(alphha - 1)*alphha*taua*exp(Aa)*exp(xia)^alphha + (Ha*alphha*exp(Aa)*exp(xia)*(Ka*exp(xia))^(alphha - 1)*(alphha - 1))/Ha^alphha - Ha^(1 - alphha)*Ka*Ka^(alphha - 2)*alphha*taua*exp(Aa)*exp(xia)^alphha*(alphha - 1)) - Lw30*(Qa*deltaT + Qa*Rka*thetta) + Lw16*exp(xia)*(delta - 1) - Ha^(1 - alphha)*Ka^(alphha - 1)*Lw39*alphha*exp(Aa)*exp(xia)^alphha + (Lw8*alphha*exp(Aa)*exp(xia)*(Ka*exp(xia))^(alphha - 1)*(alphha - 1))/(Ca^sigmma*Ha^alphha) + (Ha^(1 - alphha)*Ka^(alphha - 2)*Lw19*alphha*exp(Aa)*exp(xia)^alphha*(alphha - 1)*(taua - 1))/Qa) - Lw36*Qa - Lw16; 
    focW25_ss = Lw34*Qb*kappab - betta*(Lw11*(Qb*deltaT + Qb*Rkb*(thetta - 1) - Hb^(1 - alphha)*Kb^(alphha - 1)*alphha*taub*exp(Ab)*exp(xib)^alphha + (Hb*alphha*exp(Ab)*exp(xib)*(Kb*exp(xib))^(alphha - 1)*(alphha - 1))/Hb^alphha - Hb^(1 - alphha)*Kb*Kb^(alphha - 2)*alphha*taub*exp(Ab)*exp(xib)^alphha*(alphha - 1)) - Lw31*(Qb*deltaT + Qb*Rkb*thetta) + Lw17*exp(xib)*(delta - 1) - Hb^(1 - alphha)*Kb^(alphha - 1)*Lw40*alphha*exp(Ab)*exp(xib)^alphha + (Lw9*alphha*exp(Ab)*exp(xib)*(Kb*exp(xib))^(alphha - 1)*(alphha - 1))/(Cb^sigmma*Hb^alphha) + (Hb^(1 - alphha)*Kb^(alphha - 2)*Lw20*alphha*exp(Ab)*exp(xib)^alphha*(alphha - 1)*(taub - 1))/Qb) - Lw37*Qb - Lw17; 
    focW26_ss = Lw15 + Lw42*(nc*((zetta*(Ic/Ic - 1)^2)/2 + 1) + (Ic*nc*zetta*(Ic/Ic - 1))/Ic) - betta*(Lw12*((Ic^2*zetta)/Ic^3 + (2*Ic*zetta*(Ic/Ic - 1))/Ic^2) + (Ic^2*Lw42*nc*zetta*(Ic/Ic - 1))/Ic^2) + Lw12*((2*zetta*(Ic/Ic - 1))/Ic + (Ic*zetta)/Ic^2 + (Ic^3*Lambdac*zetta)/Ic^4 + (2*Ic^2*Lambdac*zetta*(Ic/Ic - 1))/Ic^3); 
    focW27_ss = Lw16 + Lw42*(na*((zetta*(Ia/Ia - 1)^2)/2 + 1) + (Ia*na*zetta*(Ia/Ia - 1))/Ia) - Lw10*((zetta*(Ia/Ia - 1)^2)/2 - Qa + (Ia*zetta*(Ia/Ia - 1))/Ia + 1) - betta*(Lw13*((Ia^2*zetta)/Ia^3 + (2*Ia*zetta*(Ia/Ia - 1))/Ia^2) - (Ia^2*Lw10*zetta*(Ia/Ia - 1))/Ia^2 + (Ia^2*Lw42*na*zetta*(Ia/Ia - 1))/Ia^2) + Lw13*((2*zetta*(Ia/Ia - 1))/Ia + (Ia*zetta)/Ia^2 + (Ia^3*Lambdaa*zetta)/Ia^4 + (2*Ia^2*Lambdaa*zetta*(Ia/Ia - 1))/Ia^3); 
    focW28_ss = Lw17 + Lw42*(nb*((zetta*(Ib/Ib - 1)^2)/2 + 1) + (Ib*nb*zetta*(Ib/Ib - 1))/Ib) - Lw11*((zetta*(Ib/Ib - 1)^2)/2 - Qb + (Ib*zetta*(Ib/Ib - 1))/Ib + 1) - betta*(Lw14*((Ib^2*zetta)/Ib^3 + (2*Ib*zetta*(Ib/Ib - 1))/Ib^2) - (Ib^2*Lw11*zetta*(Ib/Ib - 1))/Ib^2 + (Ib^2*Lw42*nb*zetta*(Ib/Ib - 1))/Ib^2) + Lw14*((2*zetta*(Ib/Ib - 1))/Ib + (Ib*zetta)/Ib^2 + (Ib^3*Lambdab*zetta)/Ib^4 + (2*Ib^2*Lambdab*zetta*(Ib/Ib - 1))/Ib^3); 
    focW29_ss = betta*((Lw18*(Qc*exp(xic)*(delta - 1) + Hc^(1 - alphha)*Kc^(alphha - 1)*alphha*exp(Ac)*exp(xic)^alphha*(tauc - 1)))/Qc^2 + Kc*Lw29*Rkc*thetta) - Kc*Lw35 - Lw12 + Kc*Lw29*deltaT + Kc*Lw32*kappac - (Lw18*exp(xic)*(delta - 1))/Qc; 
    focW30_ss = betta*((Lw19*(Qa*exp(xia)*(delta - 1) + Ha^(1 - alphha)*Ka^(alphha - 1)*alphha*exp(Aa)*exp(xia)^alphha*(taua - 1)))/Qa^2 - Ka*Lw10*Rka*(thetta - 1) + Ka*Lw30*Rka*thetta) - Lw13 + Lw10*(Ia - Ka*deltaT) - Ka*Lw36 + Ka*Lw30*deltaT + Ka*Lw33*kappaa - (Lw19*exp(xia)*(delta - 1))/Qa; 
    focW31_ss = betta*((Lw20*(Qb*exp(xib)*(delta - 1) + Hb^(1 - alphha)*Kb^(alphha - 1)*alphha*exp(Ab)*exp(xib)^alphha*(taub - 1)))/Qb^2 - Kb*Lw11*Rkb*(thetta - 1) + Kb*Lw31*Rkb*thetta) - Lw14 + Lw11*(Ib - Kb*deltaT) - Kb*Lw37 + Kb*Lw31*deltaT + Kb*Lw34*kappab - (Lw20*exp(xib)*(delta - 1))/Qb; 
    focW32_ss = Kc*Lw29*Qc*thetta - Lw18; 
    focW33_ss = Ka*Lw30*Qa*thetta - Ka*Lw10*Qa*(thetta - 1) - Lw19; 
    focW34_ss = Kb*Lw31*Qb*thetta - Kb*Lw11*Qb*(thetta - 1) - Lw20; 
    focW35_ss = Lw35 - Lw29 - Lw32*jc; 
    focW36_ss = Lw36 - Lw30 - Lw33*ja; 
    focW37_ss = Lw37 - Lw31 - Lw34*jb; 
    focW38_ss = Lw32*kappafa - Lw35 + Lw36*psiDa - Lw10*(psiDa + eta*(Fa*(psiDa - 1) - Fa*(psiDa - 1))*(psiDa - 1) - 1) + betta*(Lw10*(Rda*(psiDa - 1) + (thetta - 1)*(Rba - Rda + Rda*psiDa)) + Lw29*Rba*thetta - Lw30*thetta*(Rba - Rda + Rda*psiDa)) + Lw5*eta*(psiDa - 1) + Lw42*eta*na*(Fa*(psiDa - 1) - Fa*(psiDa - 1))*(psiDa - 1); 
    focW39_ss = Lw32*kappafb - Lw35 + Lw37*psiDb - Lw11*(psiDb + eta*(Fb*(psiDb - 1) - Fb*(psiDb - 1))*(psiDb - 1) - 1) + betta*(Lw11*(Rdb*(psiDb - 1) + (thetta - 1)*(Rbb - Rdb + Rdb*psiDb)) + Lw29*Rbb*thetta - Lw31*thetta*(Rbb - Rdb + Rdb*psiDb)) + Lw6*eta*(psiDb - 1) + Lw42*eta*nb*(Fb*(psiDb - 1) - Fb*(psiDb - 1))*(psiDb - 1); 
    focW40_ss = Lw35 + Lw4*eta + (Lw42*eta*nc*(2*Dc - 2*Dc))/2 - Lw29*Rdc*betta*thetta; 
    focW41_ss = Lambdac*Lw24*(jc*thetta - thetta + 1) - Fa*betta*(Lw10 - Lw10*thetta - Lw29*thetta + Lw30*thetta) - (Lambdaa*Lw22*(ja*thetta - thetta + 1))/psiDa + (Lambdaa*Lw27*(ja*thetta - thetta + 1))/psiDa; 
    focW42_ss = Lambdac*Lw25*(jc*thetta - thetta + 1) - Fb*betta*(Lw11 - Lw11*thetta - Lw29*thetta + Lw31*thetta) - (Lambdab*Lw23*(jb*thetta - thetta + 1))/psiDb + (Lambdab*Lw28*(jb*thetta - thetta + 1))/psiDb; 
    focW43_ss = Lw26*Rdc*(jc*thetta - thetta + 1) - Lw1*R - Lw4*Rdc - Lw49 + Lw24*(Rba - Rdc)*(jc*thetta - thetta + 1) + Lw25*(Rbb - Rdc)*(jc*thetta - thetta + 1) - Lw21*(Rdc - Rkc)*(jc*thetta - thetta + 1) + (Ic^2*Lw12*zetta*(Ic - Ic))/Ic^3 + (Ic^2*Lw12*betta*zetta*(Ic - Ic))/Ic^3; 
    focW44_ss = (Lw27*(ja*thetta - thetta + 1)*(Rba - Rda + Rda*psiDa))/psiDa - Lw2*R - Lw5*Rda - Lw50 - (Lw22*(ja*thetta - thetta + 1)*(Rba - Rda + Rda*psiDa - Rka*psiDa))/psiDa + (Ia^2*Lw13*zetta*(Ia - Ia))/Ia^3 + (Ia^2*Lw13*betta*zetta*(Ia - Ia))/Ia^3; 
    focW45_ss = (Lw28*(jb*thetta - thetta + 1)*(Rbb - Rdb + Rdb*psiDb))/psiDb - Lw3*R - Lw6*Rdb - Lw51 - (Lw23*(jb*thetta - thetta + 1)*(Rbb - Rdb + Rdb*psiDb - Rkb*psiDb))/psiDb + (Ib^2*Lw14*zetta*(Ib - Ib))/Ib^3 + (Ib^2*Lw14*betta*zetta*(Ib - Ib))/Ib^3; 
    focW46_ss = Lw26*(muc - 1) - Lw32*Nc; 
    focW47_ss = Lw27*(mua - 1) - Lw33*Na; 
    focW48_ss = Lw28*(mub - 1) - Lw34*Nb; 
    focW49_ss = Lw26*jc - Lw21*kappac - Lw24*kappafa - Lw25*kappafb; 
    focW50_ss = Lw27*ja - Lw22*kappaa; 
    focW51_ss = Lw28*jb - Lw23*kappab; 
    focW52_ss = -(Hc^(1 - alphha)*Kc^(alphha - 1)*Lw18*alphha*exp(Ac)*exp(xic)^alphha)/Qc; 
    focW53_ss = Ha^(1 - alphha)*Ka*Ka^(alphha - 1)*Lw10*alphha*exp(Aa)*exp(xia)^alphha - (Ha^(1 - alphha)*Ka^(alphha - 1)*Lw19*alphha*exp(Aa)*exp(xia)^alphha)/Qa + Ha^(1 - alphha)*Ka*Ka^(alphha - 1)*Lw10*alphha*betta*exp(Aa)*exp(xia)^alphha; 
    focW54_ss = Hb^(1 - alphha)*Kb*Kb^(alphha - 1)*Lw11*alphha*exp(Ab)*exp(xib)^alphha - (Hb^(1 - alphha)*Kb^(alphha - 1)*Lw20*alphha*exp(Ab)*exp(xib)^alphha)/Qb + Hb^(1 - alphha)*Kb*Kb^(alphha - 1)*Lw11*alphha*betta*exp(Ab)*exp(xib)^alphha; 

    % Auxiliar calculation
    % Compute Welf'(x) for U(x) to pin-down constant terms (in most cases constant = 0, except when the variable appears in the utility function)
    syms W Wc Wa Wb

    Wc = Cc^(1-sigmma)/(1-sigmma) - Hc^(1+pssi)/(1+pssi); 
    Wa = Ca^(1-sigmma)/(1-sigmma) - Ha^(1+pssi)/(1+pssi); 
    Wb = Cb^(1-sigmma)/(1-sigmma) - Hb^(1+pssi)/(1+pssi); 

    W = nc*Wc + na*Wa + nb*Wb;

    U_Cc = diff(W,Cc);
    U_Ca = diff(W,Ca);
    U_Cb = diff(W,Cb);
    U_Hc = diff(W,Hc);
    U_Ha = diff(W,Ha);
    U_Hb = diff(W,Hb);


    % Parameter values

        nc = 0.5;
        na = 0.25;
        nb = 0.25;
        betta = 0.99;
        thetta = 0.95;
        sigmma  = 1.02;
        pssi = 0.276;
        alphha = 0.333;
        delta = 0.025;

        zetta = 3.456; % taken from BDL16 where zetta = 1.728
        eta  = 0.0025; 
        etaD = 0.0025;
        deltaT = 0.003;
        psiDa = 1;
        psiDb = 1;

        kappac = 0.38;
        kappaa = 0.38;
        kappab = 0.38;
        kappafa = 0.38;
        kappafb = 0.38;

        rhoA = 0.85;
        rhok = 0.85;
        stdepsA = 0.007;
        stdxik = 0.005;

        phi_tau = 0.7;
        phi_cr  = 1; 

        phih = 0.85465*0;

    %Find SS for a given value of taxes = tools.     
    x = find_ss(tools);

    Cc  = x(1);  %(long_name='Consumption center')
    Ca  = x(2);  %(long_name='Consumption emergent country 1(incumbent)')
    Cb  = x(3);  %(long_name='Consumption emergent country 2(entrant)')

    Hc  = x(4);  %(long_name='(Labor) Hours worked in center')
    Ha  = x(5);  %(long_name='(Labor) Hours worked in emergent country 1')
    Hb  = x(6);  %(long_name='(Labor) Hours worked in emergent country 2')

    Yc  = x(7);  %(long_name='Output-Center')
    Ya  = x(8);  %(long_name='Output-Emergent country 1')
    Yb  = x(9);  %(long_name='Output-Emergent country 2')

    Kc  = x(10);  %(long_name='Capital-Center')
    Ka  = x(11);  %(long_name='Capital-Emergent country 1')
    Kb  = x(12);  %(long_name='Capital-Emergent country 2')

    Rkc  = x(13);  %(long_name='Gross rate of return of capital-Center')
    Rka  = x(14);  %(long_name='Gross rate of return of capital-Emergent country 1')
    Rkb  = x(15);  %(long_name='Gross rate of return of capital-Emergent country 2')

    Nc  = x(16);  %(long_name='Net worth of rep. bank-Center')
    Na  = x(17);  %(long_name='Net worth of rep. bank-Emergent country 1')
    Nb  = x(18);  %(long_name='Net worth of rep. bank-Emergent country 2')

    Fa  = x(19);  %(long_name='Foreign Lending-Emergent country 1')
    Fb  = x(20);  %(long_name='Foreign Lending-Emergent country 2')

    Dc  = x(21);  %(long_name='Deposits-Center')

    Rdc = x(22);  %(long_name='Deposits Rate-Center')
    Rda = x(23);  %(long_name='Deposits Rate-EME1')
    Rdb = x(24);  %(long_name='Deposits Rate-EME2')

    Rba  = x(25);  %(long_name='Foreign Loans rate-Emergent country 1')
    Rbb  = x(26);  %(long_name='Foreign Loans rate-Emergent country 2')

    jc  = x(27);  %(long_name='Derivative of value function of banks J-Center')
    ja  = x(28);  %(long_name='Derivative of value function of banks J-Emergent country 1')
    jb  = x(29);  %(long_name='Derivative of value function of banks J-Emergent country 2')

    muc  = x(30);  %(long_name='Lagrange multiplier of ICC-Center')
    mua  = x(31);  %(long_name='Lagrange multiplier of ICC-Emergent country 1')
    mub  = x(32);  %(long_name='Lagrange multiplier of ICC-Emergent country 2')

    Ac = 0;
    Aa = 0;
    Ab = 0;

    xic = 0;
    xia = 0;
    xib = 0;

    Qc = 1;
    Qa = 1;
    Qb = 1;

    Bc = 0;
    Ba = 0;
    Bb = 0;

    Ic = delta*Kc;
    Ia = delta*Ka;
    Ib = delta*Kb;

    R = 1/betta;

    Lambdac = betta;
    Lambdaa = betta;
    Lambdab = betta;

    %Setting value of policy tools as argument provided with function call
    if length(tools) == 3
        tauc = tools(1);%-0.2*0;
        taua = tools(2);%-0.20*0;
        taub = tools(3);%-0.20*0;
    elseif length(tools) == 2
        tauc = tools(1);
        taua = tools(2);
        taub = tools(2);
    elseif length(tools) == 1
        tauc = tools(1);
        taua = tools(1);
        taub = tools(1);
    end


    %Replace private solution (and predefined policy tools) in planner FOCs
    for i=1:54
        names{i} = strcat('focW',num2str(i),'_ss');
        eq.(names{i}) = subs(eval(names{i}));
    end

    utilArgs = ["Cc", "Ca", "Cb", "Hc", "Ha", "Hb"];
    for i = 1:length(utilArgs)
        namesAux{i} = strcat('U_',utilArgs{i}) ;
        eqAux.(namesAux{i}) = subs(eval(strcat('U_',utilArgs{i}))) ;
    end 

    %Rearrange the FOCs as A*LagrMult = -U' to solve for LagrMult vector
    temp = 'eqs = [';
    for i=1:nVars
        if i <= utilVars %Adj. of FOCs with partial utility terms (focs of x s.t. U'(x) =/= 0)
            temp = strcat(temp,' eq.(names{',num2str(i),'}) - eqAux.(namesAux{',num2str(i),'}) == - eqAux.(namesAux{',num2str(i),'}),');
        elseif i ~= nVars
            temp = strcat(temp,' eq.(names{',num2str(i),'}) == 0,');
        else
            temp = strcat(temp,' eq.(names{',num2str(i),'}) == 0];');
        end
    end

    eval(temp); %creates eqs

    [A,b] = equationsToMatrix(eqs);
    LMvars = symvar(eqs);

    A = eval(A);
    b = eval(b);

    %Calculation of estimated LagrMult of planner and estimation error
    %source of error: N Focs for N-#tools LagrMult (see L. Christiano Ramsey Handout)
    LMest = inv(A'*A)*A'*b ;
    estError = b - A*LMest ;
    errorNorm = norm(estError,1);

    % Sorting the LMs

    [~,LMvarsSort] = sort(LMvars);

    LMOut = LMest(LMvarsSort);

end 

% cell(1,51)
% for i = 1:51
% LMnames{i} = strcat('LM',num2str(i));
% end

% names(LMestOut) = LMnames;


%%%

%{
syms L26 L27 L28
test1 = - (69702817429106971*L26)/72057594037927936 - (317381044685795*L27)/35184372088832
test2 = - (17642043624503451*L27)/18014398509481984 - (5876160343671777*L28)/1125899906842624
test3 = - (17642043624503451*L28)/18014398509481984 - (5876160343671777*L26)/1125899906842624

% eqs =  [test1 == 0,
%         test2 == 1,
%         test3 == 5];

% [A,b] = equationsToMatrix(eqs)
% vars = symvar(eqs)

a = ...
[ ...
1, 2, 3, ...
]

b = [...
for i=1:3
    if i == 1
        b = [i, ...
    elseif i < 3
        i, ...
    else
        i]
    end
end


syms L26 L27 L28
test1 = - (69702817429106971*L26)/72057594037927936 - (317381044685795*L27)/35184372088832
test2 = - (17642043624503451*L27)/18014398509481984 - (5876160343671777*L28)/1125899906842624
test3 = - (17642043624503451*L28)/18014398509481984 - (5876160343671777*L26)/1125899906842624

temp = 'eqs = ['
for i=1:3
    if i < 3
        temp = strcat(temp,'test',num2str(i),' == 0',",")
    else
        temp = strcat(temp,'test',num2str(i),' == 1',"]")
    end
end

eval(temp) %creates eqs

[A,b] = equationsToMatrix(eqs)
vars = symvar(eqs)

%}

