function L = find_L(L0,psi,w_ss,betta,hh,Y_over_K,L_over_K,delta_ss,tau1_M,tau2_M,sig,chi,varphi)
L = csolve(@myfun,L0,[],1e-8,100,psi,w_ss,betta,hh,Y_over_K,L_over_K,delta_ss,tau1_M,tau2_M,sig,chi,varphi);
end
    
function L = myfun(L,psi,w_ss,betta,hh,Y_over_K,L_over_K,delta_ss,tau1_M,tau2_M,sig,chi,varphi)
 L =  (w_ss*((1-betta*hh) * ((1-hh)*((Y_over_K * (L/L_over_K))  - (delta_ss*(L/L_over_K))  -  tau1_M*(psi*L/L_over_K) + tau2_M*(psi*L/L_over_K)^2   ))^(-sig))/chi)^(1/varphi);

end

