function [y_,yhat3]=evaluation_block_3ord_NP(yhat3,u,i,order_var,ipred,ghx,ghu,ghxx,ghuu,ghxu,ghs2,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,threads)           
            [gyy, err] = A_times_B_kronecker_C(ghxx,yhat3,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            [guu, err] = A_times_B_kronecker_C(ghuu,u,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            [gyu, err] = A_times_B_kronecker_C(ghxu,yhat3,u,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            y2a = kron_d(yhat3,yhat3);
            [gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat3,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            u2a = kron_d(u,u);
            [guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            yu = kron_d(yhat3,u);
            [gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat3,yu,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            [gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
%             [gyy12, err] = A_times_B_kronecker_C(ghxx,yhat3,yhat3,threads);
%             mexErrCheck('A_times_B_kronecker_C', err);
            % NOte that not all terms appear in the non-pruned version e.g.
            % gyy12
%             yhat3 = ghx*yhat3 + 1/6*(gyyy + guuu) + 3/6*(gyyu + gyuu  ...
%                     + gyy12*0 + ghxss*yhat3 + ghuss*u)+ghu*u  + (gyy + guu + 2*gyu + ghs2)/2;
            yhat3 = ghx*yhat3 + 1/6*(gyyy + guuu) + 3/6*(gyyu + gyuu  ...
                    + ghxss*yhat3 + ghuss*u)+ghu*u  + (gyy + guu+ 2*gyu + ghs2)/2;
            y_(order_var,1) =  yhat3;
            yhat3 = yhat3(ipred);
end