function [ch,h,i,d,lambda_h,m,R,pi,ks,ls,qk,w,hs,ys,ps,div_s,kp,lp_bar,lp,hp,yp,pp,div_p,y,Rs,Rp,rb,b,Xb,Rp_bar,omegas,lambda_c,...
div_cb,div_sb,k,pr,dwl,utility,totutility, F_xi,int_xi,lls,yys,pw,bg,Rg]=solve_total(tau_p,GAMMA_SS)
parameter=calibration;
BETA=parameter(1); DELTA=parameter(2); THETA_S=parameter(3); THETA_P=parameter(4);
ALPHA=parameter(5); As=parameter(6); Ap=parameter(7); ETA=parameter(8); GAMMA=GAMMA_SS;
TAU_S=parameter(10); TAU_P=tau_p;SIGMA=parameter(12); OMEGAI=parameter(13); MU=parameter(14);
EPSILON=parameter(16);PI=parameter(17); RL=parameter(18); 
KM=parameter(19); BETA_B=parameter(20); BETA_E=parameter(21);PHI=parameter(22); 
ALPHA_S=parameter(23);OMEGA_P=parameter(25); OMEGAM=parameter(32);ALPHA_P=parameter(24);
KAPPA=parameter(26); XI=parameter(27);KAPPA_M=parameter(28);A=1;BETA_SB=parameter(29);
ALPHA_SB=parameter(30);ZETAL=parameter(31);ZETA=parameter(15);

pi=PI;
[lcr,fval,cflagr]=fsolve(@(x) solve_steady(x(1),x(2),x(3),x(4),x(5),PI,DELTA,THETA_S,THETA_P,ALPHA,As,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA_E,BETA,EPSILON,BETA_SB,ALPHA_SB,ZETAL),[7.2,6.12,0.924,0.22,0.39],optimset('TolFun',1e-10));
lp=lcr(1);ls=lcr(2);omegas=lcr(3);hp=lcr(4);pp=lcr(5);
[y,yp,ys,ch,h,d,m,ks,w,div_s,hs,kp,lp_bar,div_p,ps,Rs,Rp,Rp_bar,rb,Xb,F_xi,div_cb,div_sb,b,k,i,dwl,totutility,bg]=solve_star(lp,ls,omegas,hp,pp,PI,DELTA,THETA_S,THETA_P,ALPHA,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA,EPSILON,KAPPA_M,BETA_SB,ALPHA_SB,ZETAL);
R=1/BETA*PI;pw=(EPSILON-1)/EPSILON;lambda_c=1-BETA_B*R/PI;Rg=R;
lambda_h=1/ch;utility=log(ch)-KAPPA*h^(1+ETA)/(1+ETA)+KAPPA_M*log(m);lls=lp/ls;
yys=yp/ys;qk=1;pr=F_xi;
int_xi=KM/(KM-1)*(OMEGAM-OMEGAM^KM*omegas^(1-KM));
