function beta=solvebeta(alpha,delta,ggammabar,lambdak,thetabar,Rbar,zeta,qlLeY,KY);
syms betadm real;
assume(0<betadm<1);
mubmuedm=(ggammabar-betadm*Rbar)/(Rbar*(ggammabar-betadm*(1-zeta)));
phidm=qlLeY*(1-betadm-zeta*thetabar*ggammabar*(ggammabar-betadm*Rbar)/(Rbar*(ggammabar-betadm*(1-zeta))))/(betadm*alpha);
beta=solve(KY==betadm*alpha*(1-phidm)/(1-betadm/lambdak*(mubmuedm*zeta*thetabar*ggammabar/betadm+1-delta)));
beta=double(beta);
% options=optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);
% beta=fsolve(@(beta)beta*alpha*(1-(qlLeY*(1-beta-zeta*thetabar*ggammabar*(ggammabar-beta*Rbar)/(Rbar*(ggammabar-beta*(1-zeta))))/(beta*alpha)))/(1-beta/lambdak*(((ggammabar-beta*Rbar)/(Rbar*(ggammabar-beta*(1-zeta))))*zeta*thetabar*ggammabar/beta+1-delta))-KY,1,options)




% m=zeta*thetabar*ggammabar/(ggammabar-betadm*(1-zeta));
% u=m*ggammabar/Rbar;
%betadummy=solve((k+l)/(alpha+l*(1+m/(ggammabar-betadummy*u))+(k/lambdak)*(1+m/(ggammabar-betadummy*u)-delta))-betadummy==0);
%options=optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);
%beta=fsolve(@(beta)(k+l)/(alpha+l*(1+m/(ggammabar-beta*u))+(k/lambdak)*(1+m/(ggammabar-beta*u)-delta))-beta,1,options);