function [L1temp,L2temp,T1temp,T2temp,H1temp,H2temp,...
          Ytemp,atYtemp,Atemp,Ctemp,Vtemp,FVALtemp,CORcortemp,Utemp,LAM1temp,LAM2temp]=...
          solve_E2_1(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);
tau_gamma= TAU_GAMMA(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);

norm_switch=0;
f = @(LT1)foc_int(LT1,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
LT1_0=[L1_star_ftemp T1_star_ftemp+0.01];
if kttp1==2
    LT1_0=[0.9*L1_star_ftemp 0.9*T1_star_ftemp+0.01];
end
[LT,fval,resid,exitflag]= lsqnonlin(f,LT1_0,[0 0],[inf inf],options); 
if fval>10^-5
    if abs(resid(1))>abs(resid(2)) 
        norm_switch=1;
    else
        norm_switch=2;
    end
    f = @(LT1)foc_int(LT1,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
    [LT,fval,resid,exitflag]= lsqnonlin(f,LT1_0,[0 0],[inf inf],options); 
end
L1temp = LT(1);
T1temp = LT(2);
try
    L2temp = ((w1/w2)*(mu_l2_tild/mu_l1_tild)*L1temp^(1/psi_l1))^psi_l2;
    T2temp = ((w1/w2)*(mu_t2_tild/mu_t1_tild)*T1temp^(1/psi_t1))^psi_t2;
    H1temp = L_bar-L1temp-T1temp;
    H2temp = L_bar-L2temp-T2temp;
    Ytemp  = tau_b+w1*H1temp+w2*H2temp;
    atYtemp= chi*Ytemp^(1-taumu) - tau_gamma;
    mYtemp = chi*(1-taumu)*(Ytemp^(-taumu));
    MULtemp= (mu_l1*(L1temp^(1-1/psi_l1)) + mu_l2*(L2temp^(1-1/psi_l2)))^(-rhou_l);
    Ctemp = (mu_0*mYtemp*w1/(MULtemp*mu_l1_tild*L1temp^(-1/psi_l1)))^(eta_cp);
    Atemp = (1+r)*(Aprev + atYtemp - Ctemp); 
    Aim = abs(imag(Atemp))>0;
    EV = ppvalshort(EVfunc,Atemp*Aim+(1-Aim)*min(Atemp,A_grid(apoints)));
    Vtemp  = mu_0*(1/(1-1/eta_cp))*Ctemp^(1-1/eta_cp) ... 
                        - (1/(1-rhou_l))*(mu_l1*L1temp.^(1-1/psi_l1) + mu_l2*L2temp.^(1-1/psi_l2))^(1-rhou_l)...
                        - (1/(1-rhou_t))*(mu_t1*T1temp.^(1-1/psi_t1) + mu_t2*T2temp.^(1-1/psi_t2))^(1-rhou_t)...
                        + beta*EV;
    FVALtemp = fval;
    CORcortemp = 0;
    LAM1temp = 0;
    LAM2temp = 0;
    if fval>10^(-6)   
        H1temp = -1;
        H2temp = -1;
        Atemp = amin-1;
    end
catch
    H1temp = -1;
    H2temp = -1;
    Atemp = amin-1;
end

%% Corners start here
if H1temp<1e-9 || H2temp<1e-9 || Atemp<0
L1cor=zeros(1,7);
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));

% Calculate corner constrained: A (at amin)
if Atemp<amin
    LT1_0=[L1_star_ftemp T1_star_ftemp+0.01];
    if kttp1==2
        LT1_0=[0.9*L1_star_ftemp 0.9*T1_star_ftemp+0.01];
    end
    norm_switch=0;
    f = @(LT1)foc_corA(LT1,w1,w2,Aprev,parameters,tt,norm_switch);
    [LT,fval,resid,exitflag]= lsqnonlin(f,LT1_0,[0 0],[inf inf],options); 
    if (fval>10^-5 || abs(imag(LT(1)))>0) 
    if abs(resid(1))>abs(resid(2)) 
        norm_switch=1;
    else
        norm_switch=2;
    end
        f = @(LT1)foc_corA(LT1,w1,w2,Aprev,parameters,tt,norm_switch);
        [LT,fval,resid,exitflag]= lsqnonlin(f,LT1_0,[0 0],[inf inf],options); 
    end
    miter=1;
    LT1_0 = [0.95 0.95];
    while (fval>10^-5 || abs(imag(LT(1)))>0) && miter<11
        [LT,fval,resid,exitflag]= lsqnonlin(f,LT1_0,[0 0],[inf inf],options); 
        LT1_0 = max(LT1_0 - 0.1,0.01);
        miter=miter+1;
    end
    if abs(imag(LT(1)))>0 || LT(1)<0 || LT(2)<0
       L1cor(1)=10^-9;
       T1cor(1)=10^-9;
    else
        L1cor(1) = LT(1);
        T1cor(1) = LT(2);
    end
    if T1cor(1)>0
        try 
        L2cor(1) = ((w1/w2)*(mu_l2_tild/mu_l1_tild)*L1cor(1)^(1/psi_l1))^psi_l2;
        T2cor(1) = ((w1/w2)*(mu_t2_tild/mu_t1_tild)*T1cor(1)^(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) - tau_gamma;
            mYcor(1) = chi*(1-taumu)*(Ycor(1)^(-taumu));
            Acor(1)  = amin; 
            Ccor(1)  = Aprev + atYcor(1) - Acor(1)/(1+r);
            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)=1;
            if H1cor(1)<1e-9 || H2cor(1)<1e-9 || fval>10^-5
               Vcor(1)=-inf;
            end
        end
    catch
       Vcor(1)=-inf;
        end; 
    end
end

% Calculate corner constrained: H1
if H1temp<1e-9
    L2LAM1_0=[L2_star_ftemp 0];
    if kttp1==2
        L2LAM1_0=[0.9*L2_star_ftemp 0];
    end
    norm_switch=0;
    f = @(L2LAM1)foc_corH1(L2LAM1,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
    [L2LAM1,fval,resid,exitflag]= lsqnonlin(f,L2LAM1_0,[0 0],[inf inf],options); 
    if (fval>10^-5 || abs(imag(L2LAM1(1)))>0) 
    if abs(resid(1))>abs(resid(2)) 
        norm_switch=1;
    else
        norm_switch=2;
    end
        f = @(L2LAM1)foc_corH1(L2LAM1,w1,w2,Aprev,parameters,tt,EUfunc,norm_switch);
        [L2LAM1,fval,resid,exitflag]= lsqnonlin(f,L2LAM1_0,[0 0],[inf inf],options); 
    end
    miter=1;
    L2LAM1_0(1)=0.95;
    while (fval>10^-5 || abs(imag(L2LAM1(1)))>0) && miter<11
        [L2LAM1,fval,resid,exitflag]= lsqnonlin(f,L2LAM1_0,[0 0],[inf inf],options); 
        L2LAM1_0(1)=max(L2LAM1_0(1) - 0.1,0.01);
        miter=miter+1;
    end
    if abs(imag(L2LAM1(1)))>0
       L2cor(3)=10^-9;
       LAM1 = 0;
    else
        L2cor(3) = L2LAM1(1);
        LAM1     = L2LAM1(2);
        LAM1cor(3) = LAM1;
    end
    L1cor(3) = (((mu_l2_tild/mu_l1_tild)*(w1/w2 + LAM1))^(-psi_l1))*L2cor(3)^(psi_l1/psi_l2);
    T1cor(3) = L_bar - L1cor(3);
    if T1cor(3)>0
        try
        T2cor(3) = (((mu_t2_tild/mu_t1_tild)*(w1/w2 + LAM1))^( psi_t2))*T1cor(3)^(psi_t2/psi_t1);
        H1cor(3) = L_bar-L1cor(3)-T1cor(3);
        H2cor(3) = L_bar-L2cor(3)-T2cor(3);
        Ycor(3)  = tau_b+w1*H1cor(3)+w2*H2cor(3);
        if Ycor(3)>0
            atYcor(3)= chi*Ycor(3)^(1-taumu) - tau_gamma;
            mYcor(3) = chi*(1-taumu)*(Ycor(3)^(-taumu));
            MULtemp  = (mu_l1*(L1cor(3)^(1-1/psi_l1)) + mu_l2*(L2cor(3)^(1-1/psi_l2)))^(-rhou_l);
            Ccor(3)  = (mu_0*mYcor(3)*w2/(MULtemp*mu_l2_tild*L2cor(3)^(-1/psi_l2)))^(eta_cp);
            Acor(3)  = (1+r)*(Aprev + atYcor(3) - Ccor(3));
            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)=3;
            if H2cor(3)<1e-9 || Acor(3)<amin || fval>10^-5
               Vcor(3)=-inf;
            end
        end
        catch
            Vcor(3)=-inf;
        end
    end
end

% Calculate corner constrained: A H1
if H1temp<1e-9 && (Atemp<amin  || Acor(3)<amin)
    L2LAM1_0=[L2_star_ftemp 0];
    if kttp1==2
        L2LAM1_0=[0.9*L2_star_ftemp 0];
    end
    norm_switch=0;
    f = @(L2LAM1)foc_cor_A_H1(L2LAM1,w1,w2,Aprev,parameters,tt,norm_switch);
    [L2LAM1,fval,resid,exitflag]= lsqnonlin(f,L2LAM1_0,[0 0],[inf inf],options); 
    if (fval>10^-5 || abs(imag(L2LAM1(1)))>0) 
    if abs(resid(1))>abs(resid(2)) 
        norm_switch=1;
    else
        norm_switch=2;
    end
        f = @(L2LAM1)foc_cor_A_H1(L2LAM1,w1,w2,Aprev,parameters,tt,norm_switch);
        [L2LAM1,fval,resid,exitflag]= lsqnonlin(f,L2LAM1_0,[0 0],[inf inf],options); 
    end
    L2LAM1_0(1)=0.95;
    miter=1;
    while (fval>10^-5 || abs(imag(L2LAM1(1)))>0) && miter<11
        [L2LAM1,fval,resid,exitflag]= lsqnonlin(f,L2LAM1_0,[0 0],[inf inf],options); 
        L2LAM1_0(1)=max(L2LAM1_0(1) - 0.1,0.01);
        miter=miter+1;
    end
    if abs(imag(L2LAM1(1)))>0
       L2cor(5)=10^-9;
       LAM1 = 0;
    else
        L2cor(5) = L2LAM1(1);
        LAM1     = L2LAM1(2);
        LAM1cor(5) = LAM1;
    end
    L1cor(5) = (((mu_l2_tild/mu_l1_tild)*(w1/w2 + LAM1))^(-psi_l1))*L2cor(5)^(psi_l1/psi_l2);
    T1cor(5) = L_bar - L1cor(5);
    if T1cor(5)>0
        try
        T2cor(5) = (((mu_t2_tild/mu_t1_tild)*(w1/w2 + LAM1))^( psi_t2))*T1cor(5)^(psi_t2/psi_t1);
        H1cor(5) = L_bar-L1cor(5)-T1cor(5);
        H2cor(5) = L_bar-L2cor(5)-T2cor(5);
        Ycor(5)  = tau_b+w1*H1cor(5)+w2*H2cor(5);
        if Ycor(5)>0
            atYcor(5)= chi*Ycor(5)^(1-taumu) - tau_gamma;
            mYcor(5) = chi*(1-taumu)*(Ycor(5)^(-taumu));
            Acor(5)  = amin; 
            Ccor(5)  = Aprev + atYcor(5) - Acor(5)/(1+r);
            Aim = abs(imag(Acor(5)))>0;
            EV = ppvalshort(EVfunc,Acor(5)*Aim+(1-Aim)*min(Acor(5),A_grid(apoints)));
            Vcor(5)  = mu_0*(1/(1-1/eta_cp))*Ccor(5)^(1-1/eta_cp) ... 
                                - (1/(1-rhou_l))*(mu_l1*L1cor(5).^(1-1/psi_l1) + mu_l2*L2cor(5).^(1-1/psi_l2))^(1-rhou_l)...
                                - (1/(1-rhou_t))*(mu_t1*T1cor(5).^(1-1/psi_t1) + mu_t2*T2cor(5).^(1-1/psi_t2))^(1-rhou_t)...
                                + beta*EV;
            FVALcor(5)=fval;
            CORcor(5)=5;
            if H2cor(5)<1e-9  || fval>10^-5
               Vcor(5)=-inf;
            end
        end
        catch
            Vcor(5)=-inf;
        end
    end
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);
LAM1temp = LAM1cor(ii);
LAM2temp = LAM2cor(ii);

end

Utemp  = mu_0*Ctemp.^(-1/eta_cp);

if imag(Vtemp)>0
    L1temp = NaN; 
    T1temp = NaN;
    L2temp = NaN;
    T2temp = NaN;
    H1temp = NaN; 
    H2temp = NaN;
    Ytemp  = NaN; 
    atYtemp= NaN;
    Atemp  = NaN; 
    Ctemp  = NaN;
    FVALtemp  = NaN;
    CORcortemp  = NaN;
    LAM1temp = NaN;
    LAM2temp = NaN;
    Vtemp  = -inf;
    Utemp  = NaN;
end


end


