function [L1temp,L2temp,T1temp,T2temp,H1temp,H2temp,...
          Ytemp,atYtemp,Atemp,Ctemp,Vtemp,FVALtemp,CORcortemp,Utemp,LAM1temp,LAM2temp]=...
          solve_E2_0(Aprev,w1,w2,parameters,EUfunc,EVfunc,...
          L1_star_ftemp,T1_star_ftemp,L2_star_ftemp,A_star_ftemp,A_grid,options,tt)

eta_cp = parameters.eta_cp ;
MU_0  = parameters.mu_0  ;
psi_l1 = parameters.psi_l1 ;
psi_l2 = parameters.psi_l2 ;
psi_t1 = parameters.psi_t1 ;
psi_t2 = parameters.psi_t2 ;
rhou_l = parameters.rhou_l ;
rhou_t = parameters.rhou_t ;
MU_l1 = parameters.mu_l1 ;
MU_l2 = parameters.mu_l2 ;
MU_t1 = parameters.mu_t1 ;
MU_t2 = parameters.mu_t2 ;
L_bar = parameters.L_bar ;
KT    = parameters.KT;
CHI    = parameters.chi    ;
TAUMU  = parameters.taumu  ;
TAU_B  = parameters.tau_b  ;
TAU_GAMMA= parameters.tau_gamma;
r      = parameters.r  ;
amax   = parameters.amax   ;
amin   = parameters.amin   ;
beta   = parameters.beta  ;
apoints= parameters.apoints  ;
wpoints= parameters.wpoints;

kt   = KT(tt);
kttp1= KT(tt+1);

chi  = CHI(kt);
taumu= TAUMU(kt);
tau_b= TAU_B(kt);

mu_l1= MU_l1(kt);
mu_t1= MU_t1(kt);
mu_l2= MU_l2(kt);
mu_t2= MU_t2(kt);
mu_0 = MU_0(kt);

mu_l1_tild = mu_l1*(1/psi_l1-1);
mu_l2_tild = mu_l2*(1/psi_l2-1);
mu_t1_tild = mu_t1*(1/psi_t1-1);
mu_t2_tild = mu_t2*(1/psi_t2-1);

L1cor=zeros(1,4);
L2cor=zeros(size(L1cor));
T1cor=zeros(size(L1cor));
T2cor=zeros(size(L1cor));
H1cor=zeros(size(L1cor));
H2cor=zeros(size(L1cor));
Ycor =zeros(size(L1cor));
atYcor=zeros(size(L1cor));
mYcor =zeros(size(L1cor));
Acor =zeros(size(L1cor));
Ccor =zeros(size(L1cor));
Vcor =ones(size(L1cor))*(-inf);
FVALcor = zeros(size(L1cor));
CORcor = zeros(size(L1cor));
LAM1cor = zeros(size(L1cor));
LAM2cor = zeros(size(L1cor));

% Baseline constrained: H2
L1LAM2_0=[L1_star_ftemp 0];
if kttp1==2
    L1LAM2_0=[0.9*L1_star_ftemp 0];
end
norm_switch=0;
f = @(L1LAM2)foc_corH2(L1LAM2,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
[L1LAM2,fval,resid,exitflag]= lsqnonlin(f,L1LAM2_0,[0 -inf],[inf inf],options); 
if (fval>10^-5 || abs(imag(L1LAM2(1)))>0) 
if abs(resid(1))>abs(resid(2)) 
    norm_switch=1;
else
    norm_switch=2;
end
    f = @(L1LAM2)foc_corH2(L1LAM2,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
    [L1LAM2,fval,resid,exitflag]= lsqnonlin(f,L1LAM2_0,[0 -inf],[inf inf],options); 
end
miter=1;
L1LAM2_0(1)=0.95;
while (fval>10^-5 || abs(imag(L1LAM2(1)))>0) && miter<11
    [L1LAM2,fval,resid,exitflag]= lsqnonlin(f,L1LAM2_0,[0 -inf],[inf inf],options); 
    L1LAM2_0(1)=max(L1LAM2_0(1) - 0.1,0.01);
    miter=miter+1;
end
if abs(imag(L1LAM2(1)))>0 
   L1cor(1)=10^-9;
   LAM2 = 0;
else
    L1cor(1) = L1LAM2(1);
    LAM2     = L1LAM2(2);
    LAM2cor(1) = LAM2;
end
L2cor(1) = (((mu_l1_tild/mu_l2_tild)*(w2/w1 + LAM2))^(-psi_l2))*L1cor(1)^(psi_l2/psi_l1);
T2cor(1) = L_bar - L2cor(1);
if T2cor(1)>0
    try
    T1cor(1) = (((mu_t1_tild/mu_t2_tild)*(w2/w1 + LAM2))^( psi_t1))*T2cor(1)^(psi_t1/psi_t2);
    H1cor(1) = L_bar-L1cor(1)-T1cor(1);
    H2cor(1) = L_bar-L2cor(1)-T2cor(1);
    Ycor(1)  = tau_b+w1*H1cor(1)+w2*H2cor(1);
    if Ycor(1)>0
        atYcor(1)= chi*Ycor(1)^(1-taumu);
        mYcor(1) = chi*(1-taumu)*(Ycor(1)^(-taumu));
        MULtemp  = (mu_l1*(L1cor(1)^(1-1/psi_l1)) + mu_l2*(L2cor(1)^(1-1/psi_l2)))^(-rhou_l);
        Ccor(1)  = (mu_0*mYcor(1)*w1/(MULtemp*mu_l1_tild*L1cor(1)^(-1/psi_l1)))^(eta_cp);
        Acor(1)  = (1+r)*(Aprev + atYcor(1) - Ccor(1));
        Aim = abs(imag(Acor(1)))>0;
        EV = ppvalshort(EVfunc,Acor(1)*Aim+(1-Aim)*min(Acor(1),A_grid(apoints)));
        Vcor(1)  = mu_0*(1/(1-1/eta_cp))*Ccor(1)^(1-1/eta_cp) ... 
                            - (1/(1-rhou_l))*(mu_l1*L1cor(1).^(1-1/psi_l1) + mu_l2*L2cor(1).^(1-1/psi_l2))^(1-rhou_l)...
                            - (1/(1-rhou_t))*(mu_t1*T1cor(1).^(1-1/psi_t1) + mu_t2*T2cor(1).^(1-1/psi_t2))^(1-rhou_t)...
                            + beta*EV;
        FVALcor(1)=fval;
        CORcor(1)=11;
        if H1cor(1)<1e-9 || Acor(1)<amin || fval>10^-5
           Vcor(1)=-inf;
        end
    end
    catch
        Vcor(1)=-inf;
    end
end

% Calculate case constrained: A H2
if Acor(1)<amin || abs(imag(L1LAM2(1)))>0
    L1LAM2_0=[L1_star_ftemp 0];
    if kttp1==2
        L1LAM2_0=[0.9*L1_star_ftemp 0];
    end
    norm_switch=0;
    f = @(L1LAM2)foc_cor_A_H2(L1LAM2,w1,w2,Aprev,parameters,tt,norm_switch);
    [L1LAM2,fval,resid,exitflag]= lsqnonlin(f,L1LAM2_0,[0 -inf],[inf inf],options); 
    if (fval>10^-5 || abs(imag(L1LAM2(1)))>0) 
    if abs(resid(1))>abs(resid(2)) 
        norm_switch=1;
    else
        norm_switch=2;
    end
        f = @(L1LAM2)foc_cor_A_H2(L1LAM2,w1,w2,Aprev,parameters,tt,norm_switch);
        [L1LAM2,fval,resid,exitflag]= lsqnonlin(f,L1LAM2_0,[0 -inf],[inf inf],options); 
    end
    miter=1;
    L1LAM2_0(1)=0.95;
    while (fval>10^-5 || abs(imag(L1LAM2(1)))>0) && miter<11
        [L1LAM2,fval,resid,exitflag]= lsqnonlin(f,L1LAM2_0,[0 -inf],[inf inf],options); 
        L1LAM2_0(1)=max(L1LAM2_0(1) - 0.1,0.01);
        miter=miter+1;
    end
    if abs(imag(L1LAM2(1)))>0
       L1cor(2)=10^-9;
       LAM2 = 0;
    else
        L1cor(2) = L1LAM2(1);
        LAM2     = L1LAM2(2);
        LAM2cor(2) = LAM2;
    end
    L2cor(2) = (((mu_l1_tild/mu_l2_tild)*(w2/w1 + LAM2))^(-psi_l2))*L1cor(2)^(psi_l2/psi_l1);
    T2cor(2) = L_bar - L2cor(2);
    if T2cor(2)>0
        try
        T1cor(2) = (((mu_t1_tild/mu_t2_tild)*(w2/w1 + LAM2))^( psi_t1))*T2cor(2)^(psi_t1/psi_t2);
        H1cor(2) = L_bar-L1cor(2)-T1cor(2);
        H2cor(2) = L_bar-L2cor(2)-T2cor(2);
        Ycor(2)  = tau_b+w1*H1cor(2)+w2*H2cor(2);
        if Ycor(2)>0
            atYcor(2)= chi*Ycor(2)^(1-taumu);
            mYcor(2) = chi*(1-taumu)*(Ycor(2)^(-taumu));
            Acor(2)  = amin; 
            Ccor(2)  = Aprev + atYcor(2) - Acor(2)/(1+r);
            Aim = abs(imag(Acor(2)))>0;
            EV = ppvalshort(EVfunc,Acor(3)*Aim+(1-Aim)*min(Acor(2),A_grid(apoints)));
            Vcor(2)  = mu_0*(1/(1-1/eta_cp))*Ccor(2)^(1-1/eta_cp) ... 
                                - (1/(1-rhou_l))*(mu_l1*L1cor(2).^(1-1/psi_l1) + mu_l2*L2cor(2).^(1-1/psi_l2))^(1-rhou_l)...
                                - (1/(1-rhou_t))*(mu_t1*T1cor(2).^(1-1/psi_t1) + mu_t2*T2cor(2).^(1-1/psi_t2))^(1-rhou_t)...
                                + beta*EV;
            FVALcor(2)=fval;
            CORcor(2)=12;
            if H1cor(2)<0 || fval>10^-5
               Vcor(2)=-inf;
            end
        end
        catch
            Vcor(2)=-inf;
        end
    end
end

% Calculate case constrained: H1 H2
if Aprev>0 && H1cor(1)<1e-9 
    L1L2A_0=[L1_star_ftemp L2_star_ftemp A_star_ftemp];
    if kttp1==2
        L1L2A_0=[0.9*L1_star_ftemp 0.9*L2_star_ftemp A_star_ftemp];
    end
    norm_switch=0;
    f = @(L1L2A)foc_corH1H2(L1L2A,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
    [L1L2A,fval,resid,exitflag]= lsqnonlin(f,L1L2A_0,[0 0 -inf],[inf inf inf],options); 
    if (fval>10^-5 || abs(imag(L1L2A(1)))>0) 
        if abs(resid(1))>0.001 
            norm_switch=1;
        end
        if abs(resid(2))>0.001 
            norm_switch=2;
        end
        if abs(resid(3))>0.001 
            norm_switch=2;
        end
        if abs(resid(1))>0.001 && abs(resid(2))>0.001 
            norm_switch=12;
        end
        if abs(resid(1))>0.001 && abs(resid(3))>0.001 
            norm_switch=13;
        end
        if abs(resid(2))>0.001 && abs(resid(3))>0.001 
            norm_switch=23;
        end
        if abs(resid(1))>0.001 && abs(resid(2))>0.001 && abs(resid(3))>0.001 
            norm_switch=123;
        end
        f = @(L1L2A)foc_corH1H2(L1L2A,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
        [L1L2A,fval,resid,exitflag]= lsqnonlin(f,L1L2A_0,[0 0 -inf],[inf inf inf],options); 
    end
    if fval>10^(-6)
        L1L2A_0=[0.1 0.1 A_grid(2)];
        [L1L2A,fval,resid,exitflag]= lsqnonlin(f,L1L2A_0,[0 0 -inf],[inf inf inf],options); 
    end
    if abs(imag(L1L2A(1)))>0
        L1cor(3) = 10^-9;
        L2cor(3) = 10^-9;
        Acor(3)  = Aprev;
    else
        L1cor(3) = L1L2A(1);
        L2cor(3) = L1L2A(2);
        Acor(3)  = L1L2A(3);
    end
    T1cor(3) = L_bar - L1cor(3);
    T2cor(3) = L_bar - L2cor(3);
    if T1cor(3)>0 && T2cor(3)>0
        try
        H1cor(3) = 0;
        H2cor(3) = 0;
        Ycor(3)  = tau_b + w1*H1cor(3) + w2*H2cor(3);
        if Ycor(3)>0
            atYcor(3)= chi*Ycor(3)^(1-taumu);
            Ccor(3)  = atYcor(3) + Aprev-Acor(3)/(1+r);
            Aim = abs(imag(Acor(3)))>0;
            EV = ppvalshort(EVfunc,Acor(3)*Aim+(1-Aim)*min(Acor(3),A_grid(apoints)));
            Vcor(3)  = mu_0*(1/(1-1/eta_cp))*Ccor(3)^(1-1/eta_cp) ... 
                                - (1/(1-rhou_l))*(mu_l1*L1cor(3).^(1-1/psi_l1) + mu_l2*L2cor(3).^(1-1/psi_l2))^(1-rhou_l)...
                                - (1/(1-rhou_t))*(mu_t1*T1cor(3).^(1-1/psi_t1) + mu_t2*T2cor(3).^(1-1/psi_t2))^(1-rhou_t)...
                                + beta*EV;
            FVALcor(3)=fval;
            CORcor(3)=13;
            if Acor(3)<amin || Ccor(3)<=0 || fval>10^-5
               Vcor(3)=-inf;
            end
        end
        catch
           Vcor(3)=-inf;
        end
    end
end

% Calculate case constrained: A H1 H2
L1L2_0=[0.1 0.1];
norm_switch=0;
f = @(L1L2)foc_cor_A_H1H2(L1L2,w1,w2,Aprev,parameters,tt,norm_switch);
[L1L2,fval,resid,exitflag]= lsqnonlin(f,L1L2_0,[0 0],[inf inf],options); 
L1cor(4) = L1L2(1);
L2cor(4) = L1L2(2);
T1cor(4) = L_bar - L1cor(4);
T2cor(4) = L_bar - L2cor(4);
try
H1cor(4) = 0;
H2cor(4) = 0;
Ycor(4)  = tau_b + w1*H1cor(4) + w2*H2cor(4);
if Ycor(4)>0
    atYcor(4)= chi*Ycor(4)^(1-taumu);
    Acor(4)  = amin;
    Ccor(4)  = Aprev + atYcor(4) - Acor(4)/(1+r);
    EV       = ppvalshort(EVfunc,Acor(4));
    Vcor(4)  = mu_0*(1/(1-1/eta_cp))*Ccor(4)^(1-1/eta_cp) ... 
                        - (1/(1-rhou_l))*(mu_l1*L1cor(4).^(1-1/psi_l1) + mu_l2*L2cor(4).^(1-1/psi_l2))^(1-rhou_l)...
                        - (1/(1-rhou_t))*(mu_t1*T1cor(4).^(1-1/psi_t1) + mu_t2*T2cor(4).^(1-1/psi_t2))^(1-rhou_t)...
                        + beta*EV;
    FVALcor(4)=fval;
    CORcor(4)=14;
end
catch
    Vcor(4)=-inf;
end

%% Choose the maximum within corners
[~,ii]=max(Vcor);
L1temp = L1cor(ii); 
T1temp = T1cor(ii);
L2temp = L2cor(ii);
T2temp = T2cor(ii);
H1temp = H1cor(ii); 
H2temp = H2cor(ii);
Ytemp  = Ycor(ii); 
atYtemp= atYcor(ii);
Atemp  = Acor(ii); 
Ctemp  = Ccor(ii);
Vtemp  = Vcor(ii);
FVALtemp  = FVALcor(ii);
CORcortemp  = CORcor(ii);
Utemp  = mu_0*Ctemp.^(-1/eta_cp);
LAM1temp = LAM1cor(ii);
LAM2temp = LAM2cor(ii);

end



