function output_vf = solve_value_function_kid_lbar(parameters,wage_grids_outputs)

%% Read parameters
TT     = parameters.TT      ;
eta_cp = parameters.eta_cp ;
MU_0  = parameters.mu_0  ;
psi_l1 = parameters.psi_l1 ;
psi_l2 = parameters.psi_l2 ;
rhou_l = parameters.rhou_l ;
MU_l1 = parameters.mu_l1 ;
MU_l2 = parameters.mu_l2 ;
psi_t1 = parameters.psi_t1 ;
psi_t2 = parameters.psi_t2 ;
rhou_t = parameters.rhou_t ;
MU_t1 = parameters.mu_t1 ;
MU_t2 = parameters.mu_t2 ;
amax   = parameters.amax   ;
amin   = parameters.amin   ;
apoints= parameters.apoints;
L_bar  = parameters.L_bar;
beta   = parameters.beta   ;
r      = parameters.r      ;        
div    = parameters.div;
wpoints= parameters.wpoints;
KT     = parameters.KT ;
CHI    = parameters.chi    ;
TAUMU  = parameters.taumu  ;
TAU_B  = parameters.tau_b  ;
TAU_GAMMA= parameters.tau_gamma;
 
v1_grid    = wage_grids_outputs.v1_grid   ;
v2_grid    = wage_grids_outputs.v2_grid   ;
u1_grid    = wage_grids_outputs.u1_grid   ;
u2_grid    = wage_grids_outputs.u2_grid   ;
v1_prob    = wage_grids_outputs.v1_prob   ;
v2_prob    = wage_grids_outputs.v2_prob   ;
u1_prob    = wage_grids_outputs.u1_prob   ;
u2_prob    = wage_grids_outputs.u2_prob   ;
lw_grid_1  = wage_grids_outputs.lw_grid_1 ;
lw_grid_2  = wage_grids_outputs.lw_grid_2 ;
gamma_u    = wage_grids_outputs.gamma_u   ; 
gamma_v    = wage_grids_outputs.gamma_v   ;   
sig_u2_hat = wage_grids_outputs.sig_u2_hat;
sig_v2_hat = wage_grids_outputs.sig_v2_hat;
log_w_det_1= wage_grids_outputs.log_w_det_1;
log_w_det_2= wage_grids_outputs.log_w_det_2;

vpoints= size(v1_grid,1); 
upoints= size(u1_grid,1);

%% step 1 Set up grids, policy matrices and working matrices for optimization
% dimensions:
% ==========
% 1 assets today
% 2 permanent wage husband
% 3 permanent wage wife
% 4 transitory shock husband
% 5 transitory shock husband
% 6 time

A_grid_temp1 = linspace(0,10000,10);
A_grid_temp2 = linspace(log(10500),amax,apoints-10); 
A_grid = [A_grid_temp1,exp(A_grid_temp2)]+amin;

A_grid = A_grid/4160;

% Policy matrices
A_star  = zeros(size(A_grid,2),wpoints,wpoints,upoints,upoints,TT); % Assets
L1_star = zeros(size(A_star)); % Leisure husband
L2_star = zeros(size(A_star)); % Leisure wife
T1_star = zeros(size(A_star)); % Time with children husband
T2_star = zeros(size(A_star)); % Time with children wife
FVAL    = zeros(size(L1_star)); % Keep track of function value for diagnostics

A_star_E2_1  = zeros(size(A_star)); 
L1_star_E2_1 = zeros(size(A_star)); 
L2_star_E2_1 = zeros(size(A_star)); 
T1_star_E2_1 = zeros(size(A_star)); 
T2_star_E2_1 = zeros(size(A_star)); 
FVAL_E2_1    = zeros(size(A_star)); 

A_star_E2_0  = zeros(size(A_star)); 
L1_star_E2_0 = zeros(size(A_star)); 
L2_star_E2_0 = zeros(size(A_star)); 
T1_star_E2_0 = zeros(size(A_star)); 
T2_star_E2_0 = zeros(size(A_star)); 
FVAL_E2_0    = zeros(size(A_star)); 

% Value matrices
V_star = zeros(size(A_star));   % Value function
V_star_E2_1 = zeros(size(A_star));   
V_star_E2_0 = zeros(size(A_star));   

% Additional matrices
C_star =  zeros(size(A_star));   % The value of consumption
C_star_E2_1 =  zeros(size(A_star));   % The value of consumption when wife employed
C_star_E2_0 =  zeros(size(A_star));   % The value of consumption when wife non-employed
Y_star = zeros(size(A_star)); % The value of the before tax earnings
Y_star_E2_1 = zeros(size(A_star)); 
Y_star_E2_0 = zeros(size(A_star)); 
atY_star = zeros(size(A_star)); % The value of the after tax earnings
atY_star_E2_1 = zeros(size(A_star)); 
atY_star_E2_0 = zeros(size(A_star)); 
H1_star = zeros(size(A_star)); % Time spent in market hours for the husband
H1_star_E2_1 = zeros(size(A_star)); 
H1_star_E2_0 = zeros(size(A_star)); 
H2_star = zeros(size(A_star)); % Time spent in market hours for the husband
H2_star_E2_1 = zeros(size(A_star)); 
H2_star_E2_0 = zeros(size(A_star)); 
CORNER = zeros(size(A_star)); % Records the corner which was met
CORNER_E2_1 = zeros(size(A_star)); 
CORNER_E2_0 = zeros(size(A_star)); 
U_star = zeros(size(A_star));   % Marginal utility of consumption 
U_star_E2_1 = zeros(size(A_star));   
U_star_E2_0 = zeros(size(A_star));   
LAM1_star = zeros(size(A_star)); % The multiplier on the hours of the husband constraint
LAM1_star_E2_1 = zeros(size(A_star)); 
LAM1_star_E2_0 = zeros(size(A_star)); 
LAM2_star = zeros(size(A_star)); % The multiplier on the hours of the wife constraint
LAM2_star_E2_1 = zeros(size(A_star)); 
LAM2_star_E2_0 = zeros(size(A_star)); 

options = optimoptions('lsqnonlin','Display','off','TolX',10^-5,'TolFun',10^-12,'MaxFunEvals',5000,'MaxIter',500);

%% Maximization  
%% Last period - no choice of assets
tt = TT;
parfor i=1:apoints;
    for j=1:wpoints;
        for k=1:wpoints;
            for l=1:upoints;
                for m=1:upoints;
                    
                    %% Always use the case with no kids
                        w1 = exp(lw_grid_1(j,tt) + u1_grid(l,tt));
                        w2 = exp(lw_grid_2(k,tt) + u2_grid(m,tt));
                        Aprev=A_grid(i);
                        
                        L1_E2=zeros(1,2);
                        L2_E2=zeros(size(L1_E2));
                        T1_E2=zeros(size(L1_E2));
                        T2_E2=zeros(size(L1_E2));
                        H1_E2=zeros(size(L1_E2));
                        H2_E2=zeros(size(L1_E2));
                        Y_E2=zeros(size(L1_E2));
                        atY_E2=zeros(size(L1_E2));
                        A_E2=zeros(size(L1_E2));
                        C_E2=zeros(size(L1_E2));
                        V_E2=zeros(size(L1_E2))-inf;
                        F_E2=zeros(size(L1_E2));
                        COR_E2=zeros(size(L1_E2));
                        U_E2=zeros(size(L1_E2));
                        LAM1_E2=zeros(size(L1_E2));
                        LAM2_E2=zeros(size(L1_E2));
                        
                        % case 2: wife employed
                        parameters1 = parameters;
                        parameters1.mu_0=parameters.mu_0.*parameters.mu_0_E2;
                        [L1_E2(2),L2_E2(2),T1_E2(2),T2_E2(2),H1_E2(2),H2_E2(2),...
                        Y_E2(2),atY_E2(2),A_E2(2),C_E2(2),V_E2(2),F_E2(2),COR_E2(2),U_E2(2),...
                        LAM1_E2(2),LAM2_E2(2)]=solve_terminal_E2_1(Aprev,w1,w2,parameters1,tt);

                        % case 1: wife not employed
                        parameters1 = parameters;
                        parameters1.tau_gamma=0*parameters.tau_gamma;
                        [L1_E2(1),L2_E2(1),T1_E2(1),T2_E2(1),H1_E2(1),H2_E2(1),...
                        Y_E2(1),atY_E2(1),A_E2(1),C_E2(1),V_E2(1),F_E2(1),COR_E2(1),U_E2(1),...
                        LAM1_E2(1),LAM2_E2(1)]=solve_terminal_E2_0(Aprev,w1,w2,parameters1,tt);
                     
                        if F_E2(1)>10^-5 && F_E2(2)<10^-5   
                            V_E2(1)=-inf;
                        end
                        [~,ii]=max(V_E2);
                        
                        V_star(i,j,k,l,m,tt)= V_E2(ii);
                        L1_star(i,j,k,l,m,tt) = L1_E2(ii);
                        T1_star(i,j,k,l,m,tt) = T1_E2(ii);
                        L2_star(i,j,k,l,m,tt) = L2_E2(ii);
                        T2_star(i,j,k,l,m,tt) = T2_E2(ii);
                        A_star(i,j,k,l,m,tt)= A_E2(ii);
                        H1_star(i,j,k,l,m,tt) = H1_E2(ii);
                        H2_star(i,j,k,l,m,tt) = H2_E2(ii);
                        Y_star(i,j,k,l,m,tt)  = Y_E2(ii);
                        atY_star(i,j,k,l,m,tt)= atY_E2(ii);
                        C_star(i,j,k,l,m,tt)= C_E2(ii);
                        U_star(i,j,k,l,m,tt)= U_E2(ii);
                        FVAL(i,j,k,l,m,tt)= F_E2(ii);
                        CORNER(i,j,k,l,m,tt)= COR_E2(ii);
                        LAM1_star(i,j,k,l,m,tt)= LAM1_E2(ii);
                        LAM2_star(i,j,k,l,m,tt)= LAM2_E2(ii);

                        V_star_E2_1(i,j,k,l,m,tt)= V_E2(2);
                        L1_star_E2_1(i,j,k,l,m,tt) = L1_E2(2);
                        T1_star_E2_1(i,j,k,l,m,tt) = T1_E2(2);
                        L2_star_E2_1(i,j,k,l,m,tt) = L2_E2(2);
                        T2_star_E2_1(i,j,k,l,m,tt) = T2_E2(2);
                        A_star_E2_1(i,j,k,l,m,tt)= A_E2(2);

                        V_star_E2_0(i,j,k,l,m,tt)= V_E2(1);
                        L1_star_E2_0(i,j,k,l,m,tt) = L1_E2(1);
                        T1_star_E2_0(i,j,k,l,m,tt) = T1_E2(1);
                        L2_star_E2_0(i,j,k,l,m,tt) = L2_E2(1);
                        T2_star_E2_0(i,j,k,l,m,tt) = T2_E2(1);
                        A_star_E2_0(i,j,k,l,m,tt)= A_E2(1);
                end
            end
        end
    end
end

%% 1 to TT-1 periods
for ttt=1:TT-1
    tt = TT-ttt;
    % Assign time varying utility parameters
    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);


    %% Step1: Expectations calculations for marginal utility and for the value function
    % calculate expectations of marginal utility next period for every 
    % point on the grid of asset and F1, F2 (not u1 u2, because i.i.d).  
    % Note: dimensions are F1, F2 (exog), and A_star (endog). This will then be
    % used for the solution of the FOC
    EUtemp = zeros(apoints,wpoints,wpoints);
    EVtemp = zeros(apoints,wpoints,wpoints);

    gr1 = log_w_det_1(tt+1)-log_w_det_1(tt);
    gr2 = log_w_det_2(tt+1)-log_w_det_2(tt);
    pv = v1_prob(:,tt+1)*v2_prob(:,tt+1)';
    pu = u1_prob(:,tt+1)*u2_prob(:,tt+1)';
    grid1prime = lw_grid_1(:,tt+1);
    grid2prime = lw_grid_2(:,tt+1);

    % Prepare max and min of grids, to evaluate the off-grid F1,F2,u1,u2 at
    % the edges.
    maxtempF1 = ones(vpoints,vpoints)*max(grid1prime);
    maxtempF2 = ones(vpoints,vpoints)*max(grid2prime);
    mintempF1 = ones(vpoints,vpoints)*min(grid1prime);
    mintempF2 = ones(vpoints,vpoints)*min(grid2prime);

    maxtempu2 = ones(upoints,upoints)*max(u2_grid(:,tt+1));
    mintempu2 = ones(upoints,upoints)*min(u2_grid(:,tt+1));

    parfor k=1:apoints;
        Upar4 = squeeze(U_star(k,:,:,:,:,tt+1));
        Vpar4 = squeeze(V_star(k,:,:,:,:,tt+1));
         [EUa,EVa]=interpexp_mex(Vpar4,Upar4,wpoints,upoints,vpoints,...
         lw_grid_1,gr1,v1_grid,u1_grid,lw_grid_2,gr2,v2_grid,u2_grid,...
         maxtempF1,maxtempF2,mintempF1,mintempF2,maxtempu2,mintempu2,...
         grid1prime,grid2prime,pv,pu,gamma_v,gamma_u,tt);

        EUtemp(k,:,:) = EUa;   
        EVtemp(k,:,:) = EVa;
        
    end
    
    clear EUFUNC EVFUNC
    for j=1:wpoints;
        for k=1:wpoints;
            EUsmall = EUtemp(:,j,k);
            EVsmall = EVtemp(:,j,k);
            EUFUNC{j,k} = interp1(A_grid,EUsmall,'linear','pp');
            if EUFUNC{j,k}.coefs(1,1)>0
                EUFUNC{j,k}.coefs(1,1)=mean(EUFUNC{j,k}.coefs(2:5,1));
            end

            EVFUNC{j,k}= interp1(A_grid,EVsmall,'linear','pp');
        end
    end
    
    %% Second step: Solve for the roots 
    L1_star_f=L1_star(:,:,:,:,:,tt+1);
    T1_star_f=T1_star(:,:,:,:,:,tt+1);
    L2_star_f=L2_star(:,:,:,:,:,tt+1);
    A_star_f =A_star(:,:,:,:,:,tt+1);

    %% First solve for wife non-employed 
    parfor i=1:apoints
        parameters1 = parameters;
        parameters1.tau_gamma=0*parameters.tau_gamma;

        if kt==1
        [L1_star_E2_0_temp,T1_star_E2_0_temp,L2_star_E2_0_temp,T2_star_E2_0_temp,H1_star_E2_0_temp,H2_star_E2_0_temp,...
        Y_star_E2_0_temp,atY_star_E2_0_temp,A_star_E2_0_temp,C_star_E2_0_temp,V_star_E2_0_temp,...
        FVAL_E2_0_temp,CORNER_E2_0_temp,U_star_E2_0_temp,LAM1_star_E2_0_temp,LAM2_star_E2_0_temp, status, EU_debug]=...
        mexLab('solve_loop_E1_0', A_grid,lw_grid_1,lw_grid_2,u1_grid,u2_grid,parameters1, ... 
        tt,wpoints,upoints,i,kt,EUFUNC,EVFUNC,...
        L1_star_f,T1_star_f,L2_star_f,A_star_f);
        end

        if kt==2
        [L1_star_E2_0_temp,T1_star_E2_0_temp,L2_star_E2_0_temp,T2_star_E2_0_temp,H1_star_E2_0_temp,H2_star_E2_0_temp,...
        Y_star_E2_0_temp,atY_star_E2_0_temp,A_star_E2_0_temp,C_star_E2_0_temp,V_star_E2_0_temp,...
        FVAL_E2_0_temp,CORNER_E2_0_temp,U_star_E2_0_temp,LAM1_star_E2_0_temp,LAM2_star_E2_0_temp, status, EU_debug]=...
        mexLab('solve_loop_E1_0_nk', A_grid,lw_grid_1,lw_grid_2,u1_grid,u2_grid,parameters1, ... 
        tt,wpoints,upoints,i,kt,EUFUNC,EVFUNC,...
        L1_star_f,T1_star_f,L2_star_f,A_star_f);
        end
        
        V_star_E2_0(i,:,:,:,:,tt) =     V_star_E2_0_temp    ;
        L1_star_E2_0(i,:,:,:,:,tt) =    L1_star_E2_0_temp   ;
        T1_star_E2_0(i,:,:,:,:,tt) =    T1_star_E2_0_temp   ;
        L2_star_E2_0(i,:,:,:,:,tt) =    L2_star_E2_0_temp   ;
        T2_star_E2_0(i,:,:,:,:,tt) =    T2_star_E2_0_temp   ;
        FVAL_E2_0(i,:,:,:,:,tt) =       FVAL_E2_0_temp      ;
        A_star_E2_0(i,:,:,:,:,tt) =     A_star_E2_0_temp    ;
        C_star_E2_0(i,:,:,:,:,tt) =     C_star_E2_0_temp 	 ;
        Y_star_E2_0(i,:,:,:,:,tt) =     Y_star_E2_0_temp 	 ;
        atY_star_E2_0(i,:,:,:,:,tt) =   atY_star_E2_0_temp  ;
        H1_star_E2_0(i,:,:,:,:,tt) =    H1_star_E2_0_temp 	;
        H2_star_E2_0(i,:,:,:,:,tt) =    H2_star_E2_0_temp 	;
        CORNER_E2_0(i,:,:,:,:,tt) =     CORNER_E2_0_temp 	 ;
        U_star_E2_0(i,:,:,:,:,tt) =     U_star_E2_0_temp 	 ;
        LAM1_star_E2_0(i,:,:,:,:,tt) =  LAM1_star_E2_0_temp ;
        LAM2_star_E2_0(i,:,:,:,:,tt) =  LAM2_star_E2_0_temp ;

    end 

    %% Second solve for wife non-employed 
    parfor i=1:apoints
        if kt==1
        [L1_star_E2_1_temp,T1_star_E2_1_temp,L2_star_E2_1_temp,T2_star_E2_1_temp,H1_star_E2_1_temp,H2_star_E2_1_temp,...
        Y_star_E2_1_temp,atY_star_E2_1_temp,A_star_E2_1_temp,C_star_E2_1_temp,V_star_E2_1_temp,...
        FVAL_E2_1_temp,CORNER_E2_1_temp,U_star_E2_1_temp,LAM1_star_E2_1_temp,LAM2_star_E2_1_temp, status, EU_debug]=...
        mexLab('solve_loop_E1_1', A_grid,lw_grid_1,lw_grid_2,u1_grid,u2_grid,parameters, ... 
        tt,wpoints,upoints,i,kt,EUFUNC,EVFUNC,...
        L1_star_f,T1_star_f,L2_star_f,A_star_f);
        end
        
        if kt==2
        [L1_star_E2_1_temp,T1_star_E2_1_temp,L2_star_E2_1_temp,T2_star_E2_1_temp,H1_star_E2_1_temp,H2_star_E2_1_temp,...
        Y_star_E2_1_temp,atY_star_E2_1_temp,A_star_E2_1_temp,C_star_E2_1_temp,V_star_E2_1_temp,...
        FVAL_E2_1_temp,CORNER_E2_1_temp,U_star_E2_1_temp,LAM1_star_E2_1_temp,LAM2_star_E2_1_temp, status, EU_debug]=...
        mexLab('solve_loop_E1_1_nk', A_grid,lw_grid_1,lw_grid_2,u1_grid,u2_grid,parameters, ... 
        tt,wpoints,upoints,i,kt,EUFUNC,EVFUNC,...
        L1_star_f,T1_star_f,L2_star_f,A_star_f);
        end

        V_star_E2_1(i,:,:,:,:,tt) =     V_star_E2_1_temp    ;
        L1_star_E2_1(i,:,:,:,:,tt) =    L1_star_E2_1_temp   ;
        T1_star_E2_1(i,:,:,:,:,tt) =    T1_star_E2_1_temp   ;
        L2_star_E2_1(i,:,:,:,:,tt) =    L2_star_E2_1_temp   ;
        T2_star_E2_1(i,:,:,:,:,tt) =    T2_star_E2_1_temp   ;
        FVAL_E2_1(i,:,:,:,:,tt) =       FVAL_E2_1_temp      ;
        A_star_E2_1(i,:,:,:,:,tt) =     A_star_E2_1_temp    ;
        C_star_E2_1(i,:,:,:,:,tt) =     C_star_E2_1_temp 	 ;
        Y_star_E2_1(i,:,:,:,:,tt) =     Y_star_E2_1_temp 	 ;
        atY_star_E2_1(i,:,:,:,:,tt) =   atY_star_E2_1_temp  ;
        H1_star_E2_1(i,:,:,:,:,tt) =    H1_star_E2_1_temp 	;
        H2_star_E2_1(i,:,:,:,:,tt) =    H2_star_E2_1_temp 	;
        CORNER_E2_1(i,:,:,:,:,tt) =     CORNER_E2_1_temp 	 ;
        U_star_E2_1(i,:,:,:,:,tt) =     U_star_E2_1_temp 	 ;
        LAM1_star_E2_1(i,:,:,:,:,tt) =  LAM1_star_E2_1_temp ;
        LAM2_star_E2_1(i,:,:,:,:,tt) =  LAM2_star_E2_1_temp ;

    end

    % Finally compare the two options
    parfor i=1:apoints
    V_star_temp    = zeros(wpoints,wpoints,upoints,upoints);
    L1_star_temp   = zeros(wpoints,wpoints,upoints,upoints);
    T1_star_temp   = zeros(wpoints,wpoints,upoints,upoints);
    L2_star_temp   = zeros(wpoints,wpoints,upoints,upoints);
    T2_star_temp   = zeros(wpoints,wpoints,upoints,upoints);
    FVAL_temp      = zeros(wpoints,wpoints,upoints,upoints);
    A_star_temp    = zeros(wpoints,wpoints,upoints,upoints);
    C_star_temp 	= zeros(wpoints,wpoints,upoints,upoints);
    Y_star_temp 	= zeros(wpoints,wpoints,upoints,upoints);
    atY_star_temp  = zeros(wpoints,wpoints,upoints,upoints);
    H1_star_temp 	= zeros(wpoints,wpoints,upoints,upoints);
    H2_star_temp 	= zeros(wpoints,wpoints,upoints,upoints);
    CORNER_temp 	= zeros(wpoints,wpoints,upoints,upoints);
    U_star_temp 	= zeros(wpoints,wpoints,upoints,upoints);
    LAM1_star_temp = zeros(wpoints,wpoints,upoints,upoints);
    LAM2_star_temp = zeros(wpoints,wpoints,upoints,upoints);

    V_star_E2_0_temp     = squeeze(V_star_E2_0(i,:,:,:,:,tt));     
    L1_star_E2_0_temp    = squeeze(L1_star_E2_0(i,:,:,:,:,tt));    
    T1_star_E2_0_temp    = squeeze(T1_star_E2_0(i,:,:,:,:,tt));    
    L2_star_E2_0_temp    = squeeze(L2_star_E2_0(i,:,:,:,:,tt));    
    T2_star_E2_0_temp    = squeeze(T2_star_E2_0(i,:,:,:,:,tt));    
    FVAL_E2_0_temp       = squeeze(FVAL_E2_0(i,:,:,:,:,tt));       
    A_star_E2_0_temp     = squeeze(A_star_E2_0(i,:,:,:,:,tt));     
    C_star_E2_0_temp 	   = squeeze(C_star_E2_0(i,:,:,:,:,tt));     
    Y_star_E2_0_temp 	   = squeeze(Y_star_E2_0(i,:,:,:,:,tt));     
    atY_star_E2_0_temp   = squeeze(atY_star_E2_0(i,:,:,:,:,tt));   
    H1_star_E2_0_temp    = squeeze(H1_star_E2_0(i,:,:,:,:,tt));    
    H2_star_E2_0_temp    = squeeze(H2_star_E2_0(i,:,:,:,:,tt));    
    CORNER_E2_0_temp 	   = squeeze(CORNER_E2_0(i,:,:,:,:,tt));     
    U_star_E2_0_temp 	   = squeeze(U_star_E2_0(i,:,:,:,:,tt));     
    LAM1_star_E2_0_temp  = squeeze(LAM1_star_E2_0(i,:,:,:,:,tt));  
    LAM2_star_E2_0_temp  = squeeze(LAM2_star_E2_0(i,:,:,:,:,tt));  

    V_star_E2_1_temp     = squeeze(V_star_E2_1(i,:,:,:,:,tt));     
    L1_star_E2_1_temp    = squeeze(L1_star_E2_1(i,:,:,:,:,tt));    
    T1_star_E2_1_temp    = squeeze(T1_star_E2_1(i,:,:,:,:,tt));    
    L2_star_E2_1_temp    = squeeze(L2_star_E2_1(i,:,:,:,:,tt));    
    T2_star_E2_1_temp    = squeeze(T2_star_E2_1(i,:,:,:,:,tt));    
    FVAL_E2_1_temp       = squeeze(FVAL_E2_1(i,:,:,:,:,tt));       
    A_star_E2_1_temp     = squeeze(A_star_E2_1(i,:,:,:,:,tt));     
    C_star_E2_1_temp 	   = squeeze(C_star_E2_1(i,:,:,:,:,tt));     
    Y_star_E2_1_temp 	   = squeeze(Y_star_E2_1(i,:,:,:,:,tt));     
    atY_star_E2_1_temp   = squeeze(atY_star_E2_1(i,:,:,:,:,tt));   
    H1_star_E2_1_temp    = squeeze(H1_star_E2_1(i,:,:,:,:,tt));    
    H2_star_E2_1_temp    = squeeze(H2_star_E2_1(i,:,:,:,:,tt));    
    CORNER_E2_1_temp 	   = squeeze(CORNER_E2_1(i,:,:,:,:,tt));     
    U_star_E2_1_temp 	   = squeeze(U_star_E2_1(i,:,:,:,:,tt));     
    LAM1_star_E2_1_temp  = squeeze(LAM1_star_E2_1(i,:,:,:,:,tt));  
    LAM2_star_E2_1_temp  = squeeze(LAM2_star_E2_1(i,:,:,:,:,tt));  

    for j=1:wpoints;
        for k=1:wpoints;
            for l=1:upoints;
                for m=1:upoints;

                    if FVAL_E2_0_temp(j,k,l,m)>10^-5  
                        V_star_E2_0_temp(j,k,l,m)=-inf;
                    end
                    if FVAL_E2_1_temp(j,k,l,m)>10^-5  
                        V_star_E2_1_temp(j,k,l,m)=-inf;
                    end
                    V_star_temp(j,k,l,m)  = V_star_E2_1_temp(j,k,l,m);
                    L1_star_temp(j,k,l,m) = L1_star_E2_1_temp(j,k,l,m);
                    T1_star_temp(j,k,l,m) = T1_star_E2_1_temp(j,k,l,m);
                    L2_star_temp(j,k,l,m) = L2_star_E2_1_temp(j,k,l,m);
                    T2_star_temp(j,k,l,m) = T2_star_E2_1_temp(j,k,l,m);
                    H1_star_temp(j,k,l,m) = H1_star_E2_1_temp(j,k,l,m);
                    H2_star_temp(j,k,l,m) = H2_star_E2_1_temp(j,k,l,m);
                    Y_star_temp(j,k,l,m)  = Y_star_E2_1_temp(j,k,l,m);
                    atY_star_temp(j,k,l,m)= atY_star_E2_1_temp(j,k,l,m);
                    A_star_temp(j,k,l,m)  = A_star_E2_1_temp(j,k,l,m);
                    C_star_temp(j,k,l,m)  = C_star_E2_1_temp(j,k,l,m);
                    FVAL_temp(j,k,l,m)  = FVAL_E2_1_temp(j,k,l,m);
                    CORNER_temp(j,k,l,m)= CORNER_E2_1_temp(j,k,l,m);
                    U_star_temp(j,k,l,m)  = U_star_E2_1_temp(j,k,l,m);
                    LAM1_star_temp(j,k,l,m)  = LAM1_star_E2_1_temp(j,k,l,m);
                    LAM2_star_temp(j,k,l,m)  = LAM2_star_E2_1_temp(j,k,l,m);
                    if V_star_E2_0_temp(j,k,l,m)>V_star_E2_1_temp(j,k,l,m)
                        V_star_temp(j,k,l,m)  = V_star_E2_0_temp(j,k,l,m);
                        L1_star_temp(j,k,l,m) = L1_star_E2_0_temp(j,k,l,m);
                        T1_star_temp(j,k,l,m) = T1_star_E2_0_temp(j,k,l,m);
                        L2_star_temp(j,k,l,m) = L2_star_E2_0_temp(j,k,l,m);
                        T2_star_temp(j,k,l,m) = T2_star_E2_0_temp(j,k,l,m);
                        H1_star_temp(j,k,l,m) = H1_star_E2_0_temp(j,k,l,m);
                        H2_star_temp(j,k,l,m) = H2_star_E2_0_temp(j,k,l,m);
                        Y_star_temp(j,k,l,m)  = Y_star_E2_0_temp(j,k,l,m);
                        atY_star_temp(j,k,l,m)= atY_star_E2_0_temp(j,k,l,m);
                        A_star_temp(j,k,l,m)  = A_star_E2_0_temp(j,k,l,m);
                        C_star_temp(j,k,l,m)  = C_star_E2_0_temp(j,k,l,m);
                        FVAL_temp(j,k,l,m)    = FVAL_E2_0_temp(j,k,l,m);
                        CORNER_temp(j,k,l,m)= CORNER_E2_0_temp(j,k,l,m);
                        U_star_temp(j,k,l,m)  = U_star_E2_0_temp(j,k,l,m);
                        LAM1_star_temp(j,k,l,m)  = LAM1_star_E2_0_temp(j,k,l,m);
                        LAM2_star_temp(j,k,l,m)  = LAM2_star_E2_0_temp(j,k,l,m);
                    end
                end
            end
        end
    end    

    V_star(i,:,:,:,:,tt) =     V_star_temp   ; 
    L1_star(i,:,:,:,:,tt) =    L1_star_temp   ;
    T1_star(i,:,:,:,:,tt) =    T1_star_temp   ;
    L2_star(i,:,:,:,:,tt) =    L2_star_temp   ;
    T2_star(i,:,:,:,:,tt) =    T2_star_temp   ;
    FVAL(i,:,:,:,:,tt) =       FVAL_temp      ;
    A_star(i,:,:,:,:,tt) =     A_star_temp   ; 
    C_star(i,:,:,:,:,tt) =     C_star_temp 	; 
    Y_star(i,:,:,:,:,tt) =     Y_star_temp 	 ;
    atY_star(i,:,:,:,:,tt) =   atY_star_temp  ;
    H1_star(i,:,:,:,:,tt) =    H1_star_temp 	;
    H2_star(i,:,:,:,:,tt) =    H2_star_temp 	;
    CORNER(i,:,:,:,:,tt) =     CORNER_temp 	 ;
    U_star(i,:,:,:,:,tt) =     U_star_temp ;	 
    LAM1_star(i,:,:,:,:,tt) =  LAM1_star_temp; 
    LAM2_star(i,:,:,:,:,tt) =  LAM2_star_temp ;
    end
    
    %% Fixes few cases that GSL could not return an answer
    ind=squeeze(V_star(:,:,:,:,:,tt)<-10^30);
    [iii,jjj,kkk,lll,mmm]=ind2sub(size(ind),find(ind==1));
    if sum(ind(:))>1000
        cc=clock;
        timestamp=strcat(num2str(cc(1)),'_',num2str(cc(2)),'_',num2str(cc(3)),'_',num2str(cc(4)*100+cc(5)));
        save(strcat('debugmiss',timestamp));
        stop
    end

    % Step 1: ensure interpolation only operates on true values
    for cc=1:size(iii,1)
      ic = iii(cc); %asset
      jc = jjj(cc); %w1
      kc = kkk(cc); %w2
      lc = lll(cc); %u1
      mc = mmm(cc); %u2
      if V_star_E2_1(ic,jc,kc,lc,mc,tt)<-10^30
        V_star_E2_1(ic,jc,kc,lc,mc,tt) = -inf;
      end
      if V_star_E2_0(ic,jc,kc,lc,mc,tt)<-10^30
        V_star_E2_0(ic,jc,kc,lc,mc,tt) = -inf;
      end
    end
    
    % Step 2: Call matlab for specific points

    for cc=1:size(iii,1)
        ic = iii(cc); %asset
        jc = jjj(cc); %w1
        kc = kkk(cc); %w2
        lc = lll(cc); %u1
        mc = mmm(cc); %u2
      
        w1 = exp(lw_grid_1(jc,tt) + u1_grid(lc,tt));
        w2 = exp(lw_grid_2(kc,tt) + u2_grid(mc,tt));
        Aprev=A_grid(ic);
        EUfunc = EUFUNC{jc,kc};
        EVfunc = EVFUNC{jc,kc};

        L1_star_ftemp=L1_star_f(ic,jc,kc,lc,mc);
        T1_star_ftemp=T1_star_f(ic,jc,kc,lc,mc);
        L2_star_ftemp=L2_star_f(ic,jc,kc,lc,mc);
        A_star_ftemp=A_star_f(ic,jc,kc,lc,mc);

        if kt==1;
        parameters1 = parameters;
        parameters1.tau_gamma=0*parameters.tau_gamma;   
        [L1,L2,T1,T2,H1,H2,Y,atY,A,C,V,F,COR,U,LAM1,LAM2]=...
        solve_E2_0(Aprev,w1,w2,parameters1,EUfunc,EVfunc,...
        L1_star_ftemp,T1_star_ftemp,L2_star_ftemp,A_star_ftemp,A_grid,options,tt);
        end

        if kt==2;
        parameters1 = parameters;
        parameters1.tau_gamma=0*parameters.tau_gamma;
        [L1,L2,T1,T2,H1,H2,Y,atY,A,C,V,F,COR,U,LAM1,LAM2]=...
        solve_E2_0_nk(Aprev,w1,w2,parameters1,EUfunc,EVfunc,...
        A_grid,tt);
        end

        V_star_E2_0(ic,jc,kc,lc,mc,tt)= V;
        L1_star_E2_0(ic,jc,kc,lc,mc,tt) = L1;
        T1_star_E2_0(ic,jc,kc,lc,mc,tt) = T1;
        L2_star_E2_0(ic,jc,kc,lc,mc,tt) = L2;
        T2_star_E2_0(ic,jc,kc,lc,mc,tt) = T2;
        FVAL_E2_0(ic,jc,kc,lc,mc,tt)    = F;
        A_star_E2_0(ic,jc,kc,lc,mc,tt)  = A;
        C_star_E2_0(ic,jc,kc,lc,mc,tt) 	= C;
        Y_star_E2_0(ic,jc,kc,lc,mc,tt) 	= Y;
        atY_star_E2_0(ic,jc,kc,lc,mc,tt) = atY;	
        H1_star_E2_0(ic,jc,kc,lc,mc,tt) 	= H1;
        H2_star_E2_0(ic,jc,kc,lc,mc,tt) 	= H2;
        CORNER_E2_0(ic,jc,kc,lc,mc,tt) 	= COR;
        U_star_E2_0(ic,jc,kc,lc,mc,tt) 	= U;
        LAM1_star_E2_0(ic,jc,kc,lc,mc,tt)= LAM1;
        LAM2_star_E2_0(ic,jc,kc,lc,mc,tt)= LAM2;	

        if kt==1;
        [L1,L2,T1,T2,H1,H2,Y,atY,A,C,V,F,COR,U,LAM1,LAM2]=...
        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);
        end

        if kt==2;
        [L1,L2,T1,T2,H1,H2,Y,atY,A,C,V,F,COR,U,LAM1,LAM2]=...
        solve_E2_1_nk(Aprev,w1,w2,parameters,EUfunc,EVfunc,...
        A_grid,tt);
        end

        V_star_E2_1(ic,jc,kc,lc,mc,tt)= V;
        L1_star_E2_1(ic,jc,kc,lc,mc,tt) = L1;
        T1_star_E2_1(ic,jc,kc,lc,mc,tt) = T1;
        L2_star_E2_1(ic,jc,kc,lc,mc,tt) = L2;
        T2_star_E2_1(ic,jc,kc,lc,mc,tt) = T2;
        FVAL_E2_1(ic,jc,kc,lc,mc,tt)    = F;
        A_star_E2_1(ic,jc,kc,lc,mc,tt)  = A;
        C_star_E2_1(ic,jc,kc,lc,mc,tt) 	= C;
        Y_star_E2_1(ic,jc,kc,lc,mc,tt) 	= Y;
        atY_star_E2_1(ic,jc,kc,lc,mc,tt) = atY;	
        H1_star_E2_1(ic,jc,kc,lc,mc,tt) 	= H1;
        H2_star_E2_1(ic,jc,kc,lc,mc,tt) 	= H2;
        CORNER_E2_1(ic,jc,kc,lc,mc,tt) 	= COR;
        U_star_E2_1(ic,jc,kc,lc,mc,tt) 	= U;
        LAM1_star_E2_1(ic,jc,kc,lc,mc,tt)= LAM1;
        LAM2_star_E2_1(ic,jc,kc,lc,mc,tt)= LAM2;	
        
        if FVAL_E2_0(ic,jc,kc,lc,mc,tt)>10^-5  
            V_star_E2_0(ic,jc,kc,lc,mc,tt)=-inf;
        end
        if FVAL_E2_1(ic,jc,kc,lc,mc,tt)>10^-5  
            V_star_E2_1(ic,jc,kc,lc,mc,tt)=-inf;
        end
        V_star(ic,jc,kc,lc,mc,tt)  = V_star_E2_1(ic,jc,kc,lc,mc,tt);
        L1_star(ic,jc,kc,lc,mc,tt) = L1_star_E2_1(ic,jc,kc,lc,mc,tt);
        T1_star(ic,jc,kc,lc,mc,tt) = T1_star_E2_1(ic,jc,kc,lc,mc,tt);
        L2_star(ic,jc,kc,lc,mc,tt) = L2_star_E2_1(ic,jc,kc,lc,mc,tt);
        T2_star(ic,jc,kc,lc,mc,tt) = T2_star_E2_1(ic,jc,kc,lc,mc,tt);
        H1_star(ic,jc,kc,lc,mc,tt) = H1_star_E2_1(ic,jc,kc,lc,mc,tt);
        H2_star(ic,jc,kc,lc,mc,tt) = H2_star_E2_1(ic,jc,kc,lc,mc,tt);
        Y_star(ic,jc,kc,lc,mc,tt)  = Y_star_E2_1(ic,jc,kc,lc,mc,tt);
        atY_star(ic,jc,kc,lc,mc,tt)= atY_star_E2_1(ic,jc,kc,lc,mc,tt);
        A_star(ic,jc,kc,lc,mc,tt)  = A_star_E2_1(ic,jc,kc,lc,mc,tt);
        C_star(ic,jc,kc,lc,mc,tt)  = C_star_E2_1(ic,jc,kc,lc,mc,tt);
        FVAL(ic,jc,kc,lc,mc,tt)  = FVAL_E2_1(ic,jc,kc,lc,mc,tt);
        CORNER(ic,jc,kc,lc,mc,tt)= CORNER_E2_1(ic,jc,kc,lc,mc,tt);
        U_star(ic,jc,kc,lc,mc,tt)  = U_star_E2_1(ic,jc,kc,lc,mc,tt);
        LAM1_star(ic,jc,kc,lc,mc,tt)  = LAM1_star_E2_1(ic,jc,kc,lc,mc,tt);
        LAM2_star(ic,jc,kc,lc,mc,tt)  = LAM2_star_E2_1(ic,jc,kc,lc,mc,tt);
        if V_star_E2_0(ic,jc,kc,lc,mc,tt)>V_star_E2_1(ic,jc,kc,lc,mc,tt)
            V_star(ic,jc,kc,lc,mc,tt)  = V_star_E2_0(ic,jc,kc,lc,mc,tt);
            L1_star(ic,jc,kc,lc,mc,tt) = L1_star_E2_0(ic,jc,kc,lc,mc,tt);
            T1_star(ic,jc,kc,lc,mc,tt) = T1_star_E2_0(ic,jc,kc,lc,mc,tt);
            L2_star(ic,jc,kc,lc,mc,tt) = L2_star_E2_0(ic,jc,kc,lc,mc,tt);
            T2_star(ic,jc,kc,lc,mc,tt) = T2_star_E2_0(ic,jc,kc,lc,mc,tt);
            H1_star(ic,jc,kc,lc,mc,tt) = H1_star_E2_0(ic,jc,kc,lc,mc,tt);
            H2_star(ic,jc,kc,lc,mc,tt) = H2_star_E2_0(ic,jc,kc,lc,mc,tt);
            Y_star(ic,jc,kc,lc,mc,tt)  = Y_star_E2_0(ic,jc,kc,lc,mc,tt);
            atY_star(ic,jc,kc,lc,mc,tt)= atY_star_E2_0(ic,jc,kc,lc,mc,tt);
            A_star(ic,jc,kc,lc,mc,tt)  = A_star_E2_0(ic,jc,kc,lc,mc,tt);
            C_star(ic,jc,kc,lc,mc,tt)  = C_star_E2_0(ic,jc,kc,lc,mc,tt);
            FVAL(ic,jc,kc,lc,mc,tt)    = FVAL_E2_0(ic,jc,kc,lc,mc,tt);
            CORNER(ic,jc,kc,lc,mc,tt)= CORNER_E2_0(ic,jc,kc,lc,mc,tt);
            U_star(ic,jc,kc,lc,mc,tt)  = U_star_E2_0(ic,jc,kc,lc,mc,tt);
            LAM1_star(ic,jc,kc,lc,mc,tt)  = LAM1_star_E2_0(ic,jc,kc,lc,mc,tt);
            LAM2_star(ic,jc,kc,lc,mc,tt)  = LAM2_star_E2_0(ic,jc,kc,lc,mc,tt);
        end
        
    end
    
    % Step 3 - interpolate on still missing points
    ind=squeeze(V_star(:,:,:,:,:,tt)<-10^30);
    [iii,jjj,kkk,lll,mmm]=ind2sub(size(ind),find(ind==1));

    % V_star_E2_1 interpolations
    % Start with interpolations on the interior
    for cc=1:size(iii,1)
       ic = iii(cc); %asset
       jc = jjj(cc); %w1
       kc = kkk(cc); %w2
       lc = lll(cc); %u1
       mc = mmm(cc); %u2
       if V_star_E2_1(ic,jc,kc,lc,mc,tt)<-10^30
       if ic>1 && ic<apoints
          Ax  = [A_grid(1:ic-1),A_grid(ic+1:apoints)];
          Vy  = [V_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',V_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          L1y = [L1_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',L1_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          L2y = [L2_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',L2_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          H1y = [H1_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',H1_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          H2y = [H2_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',H2_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Fy  = [FVAL_E2_1(1:ic-1,jc,kc,lc,mc,tt)',FVAL_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Ay  = [A_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',A_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Cy  = [C_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',C_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Yy  = [Y_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',Y_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          atYy= [atY_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',atY_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Uy  = [U_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',U_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
          
          V_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Vy,A_grid(ic),'linear','extrap');
          L1_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,L1y,A_grid(ic),'linear','extrap');
          L2_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,L2y,A_grid(ic),'linear','extrap');
          H1_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,H1y,A_grid(ic),'linear','extrap');
          H2_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,H2y,A_grid(ic),'linear','extrap');
          T1_star_E2_1(ic,jc,kc,lc,mc,tt) = L_bar - L1_star_E2_1(ic,jc,kc,lc,mc,tt) - H1_star_E2_1(ic,jc,kc,lc,mc,tt);
          T2_star_E2_1(ic,jc,kc,lc,mc,tt) = L_bar - L2_star_E2_1(ic,jc,kc,lc,mc,tt) - H2_star_E2_1(ic,jc,kc,lc,mc,tt);
          FVAL_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Fy,A_grid(ic),'linear','extrap');
          A_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Ay,A_grid(ic),'linear','extrap');
          C_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Cy,A_grid(ic),'linear','extrap');
          Y_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Yy,A_grid(ic),'linear','extrap');
          atY_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,atYy,A_grid(ic),'linear','extrap');
          U_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Uy,A_grid(ic),'linear','extrap');
       end
       end
    end
    % Now can interpolate on the corners if needed
    for cc=1:size(iii,1)
       ic = iii(cc); %asset
       jc = jjj(cc); %w1
       kc = kkk(cc); %w2
       lc = lll(cc); %u1
       mc = mmm(cc); %u2

       if V_star_E2_1(ic,jc,kc,lc,mc,tt)<-10^30
       Ax  = [A_grid(1:ic-1),A_grid(ic+1:apoints)];
       Vy  = [V_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',V_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       L1y = [L1_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',L1_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       L2y = [L2_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',L2_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       H1y = [H1_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',H1_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       H2y = [H2_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',H2_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Fy  = [FVAL_E2_1(1:ic-1,jc,kc,lc,mc,tt)',FVAL_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Ay  = [A_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',A_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Cy  = [C_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',C_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Yy  = [Y_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',Y_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       atYy= [atY_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',atY_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Uy  = [U_star_E2_1(1:ic-1,jc,kc,lc,mc,tt)',U_star_E2_1(ic+1:apoints,jc,kc,lc,mc,tt)'];

       V_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Vy,A_grid(ic),'linear','extrap');
       L1_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,L1y,A_grid(ic),'linear','extrap');
       L2_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,L2y,A_grid(ic),'linear','extrap');
       H1_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,H1y,A_grid(ic),'linear','extrap');
       H2_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,H2y,A_grid(ic),'linear','extrap');
       T1_star_E2_1(ic,jc,kc,lc,mc,tt) = L_bar - L1_star_E2_1(ic,jc,kc,lc,mc,tt) - H1_star_E2_1(ic,jc,kc,lc,mc,tt);
       T2_star_E2_1(ic,jc,kc,lc,mc,tt) = L_bar - L2_star_E2_1(ic,jc,kc,lc,mc,tt) - H2_star_E2_1(ic,jc,kc,lc,mc,tt);
       FVAL_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Fy,A_grid(ic),'linear','extrap');
       A_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Ay,A_grid(ic),'linear','extrap');
       C_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Cy,A_grid(ic),'linear','extrap');
       Y_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Yy,A_grid(ic),'linear','extrap');
       atY_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,atYy,A_grid(ic),'linear','extrap');
       U_star_E2_1(ic,jc,kc,lc,mc,tt) = interp1(Ax,Uy,A_grid(ic),'linear','extrap');
       end
    end
    
    % V_star_E2_0 interpolations
    % Start with interpolations on the interior
    for cc=1:size(iii,1)
       ic = iii(cc); %asset
       jc = jjj(cc); %w1
       kc = kkk(cc); %w2
       lc = lll(cc); %u1
       mc = mmm(cc); %u2
       if V_star_E2_0(ic,jc,kc,lc,mc,tt)<-10^30
       if ic>1 && ic<apoints
          Ax  = [A_grid(1:ic-1),A_grid(ic+1:apoints)];
          Vy  = [V_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',V_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          L1y = [L1_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',L1_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          L2y = [L2_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',L2_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          H1y = [H1_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',H1_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          H2y = [H2_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',H2_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Fy  = [FVAL_E2_0(1:ic-1,jc,kc,lc,mc,tt)',FVAL_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Ay  = [A_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',A_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Cy  = [C_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',C_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Yy  = [Y_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',Y_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          atYy= [atY_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',atY_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          Uy  = [U_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',U_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
          
          V_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Vy,A_grid(ic),'linear','extrap');
          L1_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,L1y,A_grid(ic),'linear','extrap');
          L2_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,L2y,A_grid(ic),'linear','extrap');
          H1_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,H1y,A_grid(ic),'linear','extrap');
          H2_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,H2y,A_grid(ic),'linear','extrap');
          T1_star_E2_0(ic,jc,kc,lc,mc,tt) = L_bar - L1_star_E2_0(ic,jc,kc,lc,mc,tt) - H1_star_E2_0(ic,jc,kc,lc,mc,tt);
          T2_star_E2_0(ic,jc,kc,lc,mc,tt) = L_bar - L2_star_E2_0(ic,jc,kc,lc,mc,tt) - H2_star_E2_0(ic,jc,kc,lc,mc,tt);
          FVAL_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Fy,A_grid(ic),'linear','extrap');
          A_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Ay,A_grid(ic),'linear','extrap');
          C_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Cy,A_grid(ic),'linear','extrap');
          Y_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Yy,A_grid(ic),'linear','extrap');
          atY_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,atYy,A_grid(ic),'linear','extrap');
          U_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Uy,A_grid(ic),'linear','extrap');
       end
       end
    end
    % Now can interpolate on the corners if needed
    for cc=1:size(iii,1)
       ic = iii(cc); %asset
       jc = jjj(cc); %w1
       kc = kkk(cc); %w2
       lc = lll(cc); %u1
       mc = mmm(cc); %u2

       if V_star_E2_0(ic,jc,kc,lc,mc,tt)<-10^30
       Ax  = [A_grid(1:ic-1),A_grid(ic+1:apoints)];
       Vy  = [V_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',V_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       L1y = [L1_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',L1_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       L2y = [L2_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',L2_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       H1y = [H1_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',H1_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       H2y = [H2_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',H2_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Fy  = [FVAL_E2_0(1:ic-1,jc,kc,lc,mc,tt)',FVAL_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Ay  = [A_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',A_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Cy  = [C_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',C_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Yy  = [Y_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',Y_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       atYy= [atY_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',atY_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];
       Uy  = [U_star_E2_0(1:ic-1,jc,kc,lc,mc,tt)',U_star_E2_0(ic+1:apoints,jc,kc,lc,mc,tt)'];

       V_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Vy,A_grid(ic),'linear','extrap');
       L1_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,L1y,A_grid(ic),'linear','extrap');
       L2_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,L2y,A_grid(ic),'linear','extrap');
       H1_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,H1y,A_grid(ic),'linear','extrap');
       H2_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,H2y,A_grid(ic),'linear','extrap');
       T1_star_E2_0(ic,jc,kc,lc,mc,tt) = L_bar - L1_star_E2_0(ic,jc,kc,lc,mc,tt) - H1_star_E2_0(ic,jc,kc,lc,mc,tt);
       T2_star_E2_0(ic,jc,kc,lc,mc,tt) = L_bar - L2_star_E2_0(ic,jc,kc,lc,mc,tt) - H2_star_E2_0(ic,jc,kc,lc,mc,tt);
       FVAL_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Fy,A_grid(ic),'linear','extrap');
       A_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Ay,A_grid(ic),'linear','extrap');
       C_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Cy,A_grid(ic),'linear','extrap');
       Y_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Yy,A_grid(ic),'linear','extrap');
       atY_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,atYy,A_grid(ic),'linear','extrap');
       U_star_E2_0(ic,jc,kc,lc,mc,tt) = interp1(Ax,Uy,A_grid(ic),'linear','extrap');
       end
    end
    
    % Assign max of interpolated values
    for cc=1:size(iii,1)
       ic = iii(cc); %asset
       jc = jjj(cc); %w1
       kc = kkk(cc); %w2
       lc = lll(cc); %u1
       mc = mmm(cc); %u2

        if FVAL_E2_0(ic,jc,kc,lc,mc,tt)>10^-5  
            V_star_E2_0(ic,jc,kc,lc,mc,tt)=-inf;
        end
        if FVAL_E2_1(ic,jc,kc,lc,mc,tt)>10^-5  
            V_star_E2_1(ic,jc,kc,lc,mc,tt)=-inf;
        end
        V_star(ic,jc,kc,lc,mc,tt)  = V_star_E2_1(ic,jc,kc,lc,mc,tt);
        L1_star(ic,jc,kc,lc,mc,tt) = L1_star_E2_1(ic,jc,kc,lc,mc,tt);
        T1_star(ic,jc,kc,lc,mc,tt) = T1_star_E2_1(ic,jc,kc,lc,mc,tt);
        L2_star(ic,jc,kc,lc,mc,tt) = L2_star_E2_1(ic,jc,kc,lc,mc,tt);
        T2_star(ic,jc,kc,lc,mc,tt) = T2_star_E2_1(ic,jc,kc,lc,mc,tt);
        H1_star(ic,jc,kc,lc,mc,tt) = H1_star_E2_1(ic,jc,kc,lc,mc,tt);
        H2_star(ic,jc,kc,lc,mc,tt) = H2_star_E2_1(ic,jc,kc,lc,mc,tt);
        Y_star(ic,jc,kc,lc,mc,tt)  = Y_star_E2_1(ic,jc,kc,lc,mc,tt);
        atY_star(ic,jc,kc,lc,mc,tt)= atY_star_E2_1(ic,jc,kc,lc,mc,tt);
        A_star(ic,jc,kc,lc,mc,tt)  = A_star_E2_1(ic,jc,kc,lc,mc,tt);
        C_star(ic,jc,kc,lc,mc,tt)  = C_star_E2_1(ic,jc,kc,lc,mc,tt);
        FVAL(ic,jc,kc,lc,mc,tt)  = FVAL_E2_1(ic,jc,kc,lc,mc,tt);
        U_star(ic,jc,kc,lc,mc,tt)  = U_star_E2_1(ic,jc,kc,lc,mc,tt);
        LAM1_star(ic,jc,kc,lc,mc,tt)  = LAM1_star_E2_1(ic,jc,kc,lc,mc,tt);
        LAM2_star(ic,jc,kc,lc,mc,tt)  = LAM2_star_E2_1(ic,jc,kc,lc,mc,tt);
        if V_star_E2_0(ic,jc,kc,lc,mc,tt)>V_star_E2_1(ic,jc,kc,lc,mc,tt)
            V_star(ic,jc,kc,lc,mc,tt)  = V_star_E2_0(ic,jc,kc,lc,mc,tt);
            L1_star(ic,jc,kc,lc,mc,tt) = L1_star_E2_0(ic,jc,kc,lc,mc,tt);
            T1_star(ic,jc,kc,lc,mc,tt) = T1_star_E2_0(ic,jc,kc,lc,mc,tt);
            L2_star(ic,jc,kc,lc,mc,tt) = L2_star_E2_0(ic,jc,kc,lc,mc,tt);
            T2_star(ic,jc,kc,lc,mc,tt) = T2_star_E2_0(ic,jc,kc,lc,mc,tt);
            H1_star(ic,jc,kc,lc,mc,tt) = H1_star_E2_0(ic,jc,kc,lc,mc,tt);
            H2_star(ic,jc,kc,lc,mc,tt) = H2_star_E2_0(ic,jc,kc,lc,mc,tt);
            Y_star(ic,jc,kc,lc,mc,tt)  = Y_star_E2_0(ic,jc,kc,lc,mc,tt);
            atY_star(ic,jc,kc,lc,mc,tt)= atY_star_E2_0(ic,jc,kc,lc,mc,tt);
            A_star(ic,jc,kc,lc,mc,tt)  = A_star_E2_0(ic,jc,kc,lc,mc,tt);
            C_star(ic,jc,kc,lc,mc,tt)  = C_star_E2_0(ic,jc,kc,lc,mc,tt);
            FVAL(ic,jc,kc,lc,mc,tt)    = FVAL_E2_0(ic,jc,kc,lc,mc,tt);
            CORNER(ic,jc,kc,lc,mc,tt)= CORNER_E2_0(ic,jc,kc,lc,mc,tt);
            U_star(ic,jc,kc,lc,mc,tt)  = U_star_E2_0(ic,jc,kc,lc,mc,tt);
            LAM1_star(ic,jc,kc,lc,mc,tt)  = LAM1_star_E2_0(ic,jc,kc,lc,mc,tt);
            LAM2_star(ic,jc,kc,lc,mc,tt)  = LAM2_star_E2_0(ic,jc,kc,lc,mc,tt);
        end
    end
end

output_vf.A_grid   = A_grid  ;
output_vf.V_star   = V_star  ;
output_vf.U_star   = U_star  ;
output_vf.A_star = A_star  ;
output_vf.L1_star= L1_star  ;
output_vf.L2_star= L2_star  ;
output_vf.T1_star= T1_star  ;
output_vf.T2_star= T2_star  ;

output_vf.V_star_E2_1 = V_star_E2_1  ;
output_vf.A_star_E2_1 = A_star_E2_1  ;
output_vf.L1_star_E2_1= L1_star_E2_1  ;
output_vf.L2_star_E2_1= L2_star_E2_1  ;
output_vf.T1_star_E2_1= T1_star_E2_1  ;
output_vf.T2_star_E2_1= T2_star_E2_1  ;

output_vf.V_star_E2_0 = V_star_E2_0;
output_vf.A_star_E2_0 = A_star_E2_0;
output_vf.L1_star_E2_0= L1_star_E2_0;
output_vf.L2_star_E2_0= L2_star_E2_0;
output_vf.T1_star_E2_0= T1_star_E2_0;
output_vf.T2_star_E2_0= T2_star_E2_0;

%% Can be used for diagnostics
output_vf.H1_star= H1_star  ;
output_vf.H2_star= H2_star  ;
output_vf.Y_star   = Y_star  ;
output_vf.atY_star = atY_star  ;
output_vf.C_star   = C_star  ;
output_vf.FVAL     = FVAL ;
output_vf.CORNER   = CORNER ;
output_vf.LAM1_star   = LAM1_star;
output_vf.LAM2_star   = LAM2_star;
