function [EUa,EVa]=interpexp(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)

EUa = zeros(wpoints,wpoints);
EVa = zeros(wpoints,wpoints);
 
for i=1:wpoints;
    for j=1:wpoints;
        eU_u=zeros(upoints,upoints);
        eV_u=zeros(upoints,upoints);
        for uu1 = 1:upoints
            for uu2 = 1:upoints
                u4interp = squeeze(Upar4(:,:,uu1,uu2));
                v4interp = squeeze(Vpar4(:,:,uu1,uu2));
                F1prime = (lw_grid_1(i,tt) + gr1 + v1_grid(:,tt+1))*ones(1,vpoints);    
                F1prime(F1prime>maxtempF1) = maxtempF1(F1prime>maxtempF1);
                F1prime(F1prime<mintempF1) = mintempF1(F1prime<mintempF1);

                F2prime = lw_grid_2(j,tt) + gr2 + gamma_v*v1_grid(:,tt+1)*ones(1,vpoints) + ones(1,vpoints)'*v2_grid(:,tt+1)';
                F2prime(F2prime>maxtempF2) = maxtempF2(F2prime>maxtempF2);
                F2prime(F2prime<mintempF2) = mintempF2(F2prime<mintempF2);

                Uv4exp   = interp2(grid1prime,grid2prime,u4interp,F1prime,F2prime);
                Vv4exp   = interp2(grid1prime,grid2prime,v4interp,F1prime,F2prime);
                eU_u(uu1,uu2) = sum(sum(pv.*Uv4exp));
                eV_u(uu1,uu2) = sum(sum(pv.*Vv4exp));
            end
        end
        u1shock = u1_grid(:,tt+1)*ones(1,upoints);
        u2shock = gamma_u*u1shock + ones(1,upoints)'*u2_grid(:,tt+1)';
        u2shock(u2shock>maxtempu2) = maxtempu2(u2shock>maxtempu2);  
        u2shock(u2shock<mintempu2) = mintempu2(u2shock<mintempu2);

        Uu4exp = interp2(u1_grid(:,tt+1),u2_grid(:,tt+1),eU_u,u1shock,u2shock);
        Vu4exp = interp2(u1_grid(:,tt+1),u2_grid(:,tt+1),eV_u,u1shock,u2shock);

        EUa(i,j) = sum(sum(pu.*Uu4exp));   
        EVa(i,j) = sum(sum(pu.*Vu4exp));
    end
end 
 
 
end