function y = findL(LAMBDA,state0,MU)
%y = findL(LAMBDA,state0,MU) evaluates the implementability constraint of the Ramsey problem. A solution to this costraint is reached when y=0. This program is used by L.m, which searches for a value of LAMBDA that makes y=0.
%Inputs: LAMBDA (the lagrange multiplier of the implementability constraint in the Ramsey problem), state0 (the initial state of the economy), and MU (the markup).  
%Output: y, the value taken by the implementability constraint (0 is a solution)
%Calls: param.m, findvh.m
%(c) Stephanie Schmitt-Grohe and Martin Uribe

LAMBDA = abs(LAMBDA);

[BETA,THETA,D,ETA,B,A,gl,gh,phig,zl,zh,phiz,ALFA,TM]=param(MU);

% STATES
%(gh,zh)=1
%(gh,zl)=2;
%(gl,zh)=3;
%(gl,zl)=4;

g=[gh;gh;gl;gl];
z=[zh;zl;zh;zl];

g0 = g(state0);
z0 = z(state0);

OPTIONS = OPTIMSET('tolx',1e-15,'tolfun',1e-15);

Vs = sqrt(B/A);
Vmax = sqrt((1+B)/A);
P0=0; %t>0$
for j=1:4
vh = fsolve('findvh',[2.7 .2],OPTIONS,LAMBDA,g(j),z(j),MU,P0);
v(j) = Vs + (Vmax-Vs) * abs(vh(1)) / (1+abs(vh(1))); %Vs + abs(vh(1)-Vs); 
h(j) = abs(vh(2)) / (1+abs(vh(2)));
s(j)=A*v(j)+B/v(j)-2*sqrt(A*B);
c(j) = (z(j)*h(j)-g(j)) / (1+ALFA*s(j));
uc(j)=1/c(j);
uh(j)=-THETA/(1-h(j));
sp(j)=A-B/v(j)^2;
R(j)=1/(1-v(j)^2*sp(j)); %nominal interest rate
gama(j)=1+s(j)+v(j)*sp(j);
phi(j) = (1+(1-1/R(j))/v(j))/gama(j) + ALFA*s(j)/gama(j);
w(j) = (1+ETA)/ETA * z(j); %wage rate
tau(j) = 1 + uh(j)/uc(j) * gama(j) / w(j); %labor tax rate
ps(j) = tau(j)*w(j)*h(j)-g(j); %primary surplus
lhs(j,1) = uc(j)*c(j)*phi(j) + uh(j)*h(j) + uc(j)*z(j)*h(j)/gama(j)/ETA;
end

P0=1; %t=0$
vh0 = fsolve('findvh',[2.7 .2],OPTIONS,LAMBDA,g0,z0,MU,P0);

v0 = Vs + (Vmax-Vs) * abs(vh0(1)) / (1+abs(vh0(1))); %Vs + abs(vh(1)-Vs);
h0 = abs(vh0(2)) / (1+abs(vh0(2)));
s0=A*v0+B/v0-2*sqrt(A*B);
c0 = (z0*h0-g0) / (1+ALFA*s0);
uc0=1/c0;
ucc0=-1/c0^2;
uch0=0;
uhc0=0;
uh0=-THETA/(1-h0);
uhh0 = -THETA/(1-h0)^2;
sp0=A-B/v0^2;
spp0=2*B/v0^3;
R0=1/(1-v0^2*sp0);
gama0=1+s0+v0*sp0;
gamap0=2*A;
phi0 = (1+(1-1/R0)/v0)/gama0 + ALFA*s0/gama0;

lhs0 = uc0*c0*phi0 + uh0*h0 + uc0*z0*h0/gama0/ETA;

E0LHS = lhs0 + BETA * TM(state0,:) * inv(eye(size(TM))-BETA*TM) * lhs;

RHS=uc0/gama0*D;

y = E0LHS-RHS; 

y