function [y_,yhat1,yhat2,yhat3]=evaluation_block_3ord(yhat1,yhat2,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,yhat1,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,yhat1,u,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            y2a = kron_d(yhat1,yhat1);
            [gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1,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(yhat1,u);
            [gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,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,yhat1,yhat2,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
             [gyu12, err] = A_times_B_kronecker_C(ghxu,yhat2,u,threads);
            mexErrCheck('A_times_B_kronecker_C', err);
            yhat3 = ghx*yhat3 + 1/6*(gyyy + guuu) + 1/2*(gyyu + gyuu  ...
                     + ghxss*yhat1 + ghuss*u)+ gyy12+gyu12;
            yhat2 = ghx*yhat2 + 1/2*gyy +1/2* guu + gyu + 1/2*ghs2;
            yhat1 = ghx*yhat1 + ghu*u;
            y_(order_var,1) = yhat1 +yhat2 + yhat3;
            yhat1 = yhat1(ipred);
            yhat2 = yhat2(ipred);
            yhat3 = yhat3(ipred);
end