function output_vf = solve_value_function_kid_lbar_short(parameters,wage_grids_outputs,output_vf_input)

%% 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 if wife employed 
C_star_E2_0 =  zeros(size(A_star));   % The value of consumption if 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  

%% Periods with identical VF
for ttt=0:27
    tt = TT-ttt;
    V_star(:,:,:,:,:,tt) =     output_vf_input.V_star(:,:,:,:,:,tt); 
    L1_star(:,:,:,:,:,tt) =    output_vf_input.L1_star(:,:,:,:,:,tt);
    T1_star(:,:,:,:,:,tt) =    output_vf_input.T1_star(:,:,:,:,:,tt);
    L2_star(:,:,:,:,:,tt) =    output_vf_input.L2_star(:,:,:,:,:,tt);
    T2_star(:,:,:,:,:,tt) =    output_vf_input.T2_star(:,:,:,:,:,tt);
    A_star(:,:,:,:,:,tt) =     output_vf_input.A_star(:,:,:,:,:,tt); 
    C_star(:,:,:,:,:,tt) =     output_vf_input.C_star(:,:,:,:,:,tt); 
    Y_star(:,:,:,:,:,tt) =     output_vf_input.Y_star(:,:,:,:,:,tt);
    atY_star(:,:,:,:,:,tt) =   output_vf_input.atY_star(:,:,:,:,:,tt);
    H1_star(:,:,:,:,:,tt) =    output_vf_input.H1_star(:,:,:,:,:,tt);
    H2_star(:,:,:,:,:,tt) =    output_vf_input.H2_star(:,:,:,:,:,tt);
    FVAL(:,:,:,:,:,tt) =       output_vf_input.FVAL(:,:,:,:,:,tt);
    CORNER(:,:,:,:,:,tt) =     output_vf_input.CORNER(:,:,:,:,:,tt);
    U_star(:,:,:,:,:,tt) =     output_vf_input.U_star(:,:,:,:,:,tt);	 
    LAM1_star(:,:,:,:,:,tt) =  output_vf_input.LAM1_star(:,:,:,:,:,tt); 
    LAM2_star(:,:,:,:,:,tt) =  output_vf_input.LAM2_star(:,:,:,:,:,tt);

    V_star_E2_1(:,:,:,:,:,tt) =     output_vf_input.V_star_E2_1(:,:,:,:,:,tt);
    L1_star_E2_1(:,:,:,:,:,tt) =    output_vf_input.L1_star_E2_1(:,:,:,:,:,tt);
    T1_star_E2_1(:,:,:,:,:,tt) =    output_vf_input.T1_star_E2_1(:,:,:,:,:,tt);
    L2_star_E2_1(:,:,:,:,:,tt) =    output_vf_input.L2_star_E2_1(:,:,:,:,:,tt);
    T2_star_E2_1(:,:,:,:,:,tt) =    output_vf_input.T2_star_E2_1(:,:,:,:,:,tt);
    A_star_E2_1(:,:,:,:,:,tt) =     output_vf_input.A_star_E2_1(:,:,:,:,:,tt);

    V_star_E2_0(:,:,:,:,:,tt) =     output_vf_input.V_star_E2_0(:,:,:,:,:,tt);
    L1_star_E2_0(:,:,:,:,:,tt) =    output_vf_input.L1_star_E2_0(:,:,:,:,:,tt);
    T1_star_E2_0(:,:,:,:,:,tt) =    output_vf_input.T1_star_E2_0(:,:,:,:,:,tt);
    L2_star_E2_0(:,:,:,:,:,tt) =    output_vf_input.L2_star_E2_0(:,:,:,:,:,tt);
    T2_star_E2_0(:,:,:,:,:,tt) =    output_vf_input.T2_star_E2_0(:,:,:,:,:,tt);
    A_star_E2_0(:,:,:,:,:,tt) =     output_vf_input.A_star_E2_0(:,:,:,:,:,tt);

end

for ttt=28: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
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;

end