function h4d = diff4(h4,N,a1)
%DIFF4 Difference of triangular 4th order kernel in vector form
%
% H4D = DIFF4(H4,N,A) returns the triangular kernel associated to
% the generalized transfer function 
% H4D(z1,z2,z3,z4)=(1-A/z1)(1-A/z2)(1-A/z3)(1-A/z4)H4(z1,z2,z3,z4)
%
% N is the memory length (in samples) of H4 and H4D, which are vectorized

%   CBVK Toolbox
%   JHMG, PMSB 2016

h4d = zeros(size(h4));

a2 = a1*a1;
a3 = a1*a2;
a4 = a1*a3;

ind = 0;

fac1 = 1/6;
fac2 = 1/12;
fac3 = 1/24;
fac4 = 1/4;

NA2 = N^2;
NA3 = N^3;

c2N = 2*N;
c3N = 3*N;
c5N = 5*N;
c10N = 10*N;
c13N = 13*N;
c17N = 17*N;

c13Nd3 = (c13N)/3;
c5Nd2 = (c5N)/2;
c17Nd6 = (c17N)/6;

ind = ind+1;
% (1,1,1,1)
h4d(ind) = h4(ind);

ind = ind+1;
% (1,1,1,2)
h4d(ind) =  - a1*h4( ind + round( -1 ) ) + (fac4)*h4(ind);
h4d(ind) = 4*h4d(ind);

for k4=1+2:N
    ind = ind+1;
    % (1,1,1,k4)
    h4d(ind) =  - a1*(fac4)*h4( ind + round( -1 ) ) + (fac4)*h4(ind);
    h4d(ind) = 4*h4d(ind);
    
end

ind = ind+1;
% (1,1,2,2)
h4d(ind) = a2*h4( ind + round( -N ) ) - a1*2*(fac4)*h4( ind + round( 1 - N ) ) + (fac1)*h4(ind);
h4d(ind) = 6*h4d(ind);

ind = ind+1;
% (1,1,2,3)
h4d(ind) = a2*(fac4)*h4( ind + round( -N ) ) - a1*(fac1)*h4( ind + round( -1 ) ) - a1*(fac4)*h4( ind + round( 1 - N ) ) + (fac2)*h4(ind);
h4d(ind) = 12*h4d(ind);

for k4=2+2:N
    ind = ind+1;
    % (1,1,2,k4)
    h4d(ind) = a2*(fac4)*h4( ind + round( -N ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a1*(fac4)*h4( ind + round( 1 - N ) ) + (fac2)*h4(ind);
    h4d(ind) = 12*h4d(ind);
    
end

for k3=1+2:N
    ind = ind+1;
    % (1,1,k3,k3)
    h4d(ind) = a2*(fac1)*h4( ind + round( k3 - N - 2 ) ) - a1*2*(fac2)*h4( ind + round( k3 - N - 1 ) ) + (fac1)*h4(ind);
    h4d(ind) = 6*h4d(ind);
    
    
    if (k3 <= N-1)
        % (1,1,k3,k3+1)
        ind = ind+1;
        h4d(ind) = a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a1*(fac1)*h4( ind + round( -1 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    for k4=k3+0+2:N
        ind = ind+1;
        % (1,1,k3,k4)
        h4d(ind) = a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
end

ind = ind+1;
% (1,2,2,2)
h4d(ind) =  - a3*h4( ind + round( -(N*(N + 1))/2 ) ) + a2*3*(fac4)*h4( ind + round( -((N - 1)*(N + 2))/2 ) ) - a1*3*(fac1)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac4)*h4(ind);
h4d(ind) = 4*h4d(ind);

ind = ind+1;
% (1,2,2,3)
h4d(ind) =  - a3*(fac4)*h4( ind + round( -(N*(N + 1))/2 ) ) + a2*2*(fac1)*h4( ind + round( N/2 - NA2/2 - 1 ) ) - a1*(fac4)*h4( ind + round( -1 ) ) + a2*(fac4)*h4( ind + round( -((N - 1)*(N + 2))/2 ) ) - a1*2*(fac2)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac2)*h4(ind);
h4d(ind) = 12*h4d(ind);

for k4=2+2:N
    ind = ind+1;
    % (1,2,2,k4)
    h4d(ind) =  - a3*(fac4)*h4( ind + round( -(N*(N + 1))/2 ) ) + a2*2*(fac2)*h4( ind + round( N/2 - NA2/2 - 1 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) + a2*(fac4)*h4( ind + round( -((N - 1)*(N + 2))/2 ) ) - a1*2*(fac2)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac2)*h4(ind);
    h4d(ind) = 12*h4d(ind);
    
end

ind = ind+1;
% (1,2,3,3)
h4d(ind) =  - a3*(fac1)*h4( ind + round( -((N - 1)*(N + 2))/2 ) ) + a2*(fac4)*h4( ind + round( 1 - N ) ) + a2*2*(fac2)*h4( ind + round( 2 - NA2/2 - N/2 ) ) - a1*2*(fac2)*h4( ind + round( 2 - N ) ) - a1*(fac1)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac2)*h4(ind);
h4d(ind) = 12*h4d(ind);

ind = ind+1;
% (1,2,3,4)
h4d(ind) =  - a3*(fac2)*h4( ind + round( -((N - 1)*(N + 2))/2 ) ) + a2*(fac2)*h4( ind + round( 1 - N ) ) + a2*(fac1)*h4( ind + round( N/2 - NA2/2 - 1 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) + a2*(fac2)*h4( ind + round( 2 - NA2/2 - N/2 ) ) - a1*(fac2)*h4( ind + round( 2 - N ) ) - a1*(fac2)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac3)*h4(ind);
h4d(ind) = 24*h4d(ind);

for k4=3+2:N
    ind = ind+1;
    % (1,2,3,k4)
    h4d(ind) =  - a3*(fac2)*h4( ind + round( -((N - 1)*(N + 2))/2 ) ) + a2*(fac2)*h4( ind + round( 1 - N ) ) + a2*(fac2)*h4( ind + round( N/2 - NA2/2 - 1 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) + a2*(fac2)*h4( ind + round( 2 - NA2/2 - N/2 ) ) - a1*(fac2)*h4( ind + round( 2 - N ) ) - a1*(fac2)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac3)*h4(ind);
    h4d(ind) = 24*h4d(ind);
    
end

for k3=2+2:N
    ind = ind+1;
    % (1,2,k3,k3)
    h4d(ind) =  - a3*(fac1)*h4( ind + round( k3 - N/2 - NA2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) + a2*2*(fac2)*h4( ind + round( k3 - N/2 - NA2/2 - 1 ) ) - a1*2*(fac3)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac1)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac2)*h4(ind);
    h4d(ind) = 12*h4d(ind);
    
    
    if (k3 <= N-1)
        % (1,2,k3,k3+1)
        ind = ind+1;
        h4d(ind) =  - a3*(fac2)*h4( ind + round( k3 - N/2 - NA2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) + a2*(fac1)*h4( ind + round( N/2 - NA2/2 - 1 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) + a2*(fac2)*h4( ind + round( k3 - N/2 - NA2/2 - 1 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac2)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac3)*h4(ind);
        h4d(ind) = 24*h4d(ind);
        
    end
    for k4=k3+0+2:N
        ind = ind+1;
        % (1,2,k3,k4)
        h4d(ind) =  - a3*(fac2)*h4( ind + round( k3 - N/2 - NA2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) + a2*(fac2)*h4( ind + round( N/2 - NA2/2 - 1 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) + a2*(fac2)*h4( ind + round( k3 - N/2 - NA2/2 - 1 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac2)*h4( ind + round( -(N*(N - 1))/2 ) ) + (fac3)*h4(ind);
        h4d(ind) = 24*h4d(ind);
        
    end
end

for k2=1+2:N
    ind = ind+1;
    k3 = k2;
    % (1,k2,k2,k2)
    h4d(ind) =  - a3*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*3*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac4)*h4(ind);
    h4d(ind) = 4*h4d(ind);
    
    
    if (k2 <= N-1)
        % (1,k2,k2,k2+1)
        ind = ind+1;
        h4d(ind) =  - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac4)*h4( ind + round( -1 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    for k4=k2+0+2:N
        ind = ind+1;
        % (1,k2,k2,k4)
        h4d(ind) =  - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    
    if (k2 <= N-1 && k2 <= N-1)
        k3 = k2+1;
        % (1,k2,k2+1,k2+1)
        ind = ind+1;
        h4d(ind) =  - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*(fac4)*h4( ind + round( k3 - N - 2 ) ) + a2*2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*2*(fac2)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    
    if (k2 <= N-1 && k2 <= N-2)
        % (1,k2,k2+1,k2+2)
        ind = ind+1;
        h4d(ind) =  - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac3)*h4(ind);
        h4d(ind) = 24*h4d(ind);
        
    end
    for k4=k2+1+2:N
        
        if (k2 <= N-1)
            % (1,k2,k2+1,k4)
            ind = ind+1;
            h4d(ind) =  - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac3)*h4(ind);
            h4d(ind) = 24*h4d(ind);
            
        end
    end
    for k3=k2+0+2:N
        ind = ind+1;
        % (1,k2,k3,k3)
        h4d(ind) =  - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) + a2*2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
        
        if (k3 <= N-1)
            % (1,k2,k3,k3+1)
            ind = ind+1;
            h4d(ind) =  - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac3)*h4(ind);
            h4d(ind) = 24*h4d(ind);
            
        end
        for k4=k3+0+2:N
            ind = ind+1;
            % (1,k2,k3,k4)
            h4d(ind) =  - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) + (fac3)*h4(ind);
            h4d(ind) = 24*h4d(ind);
            
        end
    end
end

for k1=2:N
    ind = ind+1;
    k2 = k1;
    k3 = k1;
    % (k1,k1,k1,k1)
    h4d(ind) = a4*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*4*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*6*(fac1)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*4*(fac4)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + h4(ind);
    
    
    if (k1 <= N-1)
        % (k1,k1,k1,k1+1)
        ind = ind+1;
        h4d(ind) = a4*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*3*(fac1)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*3*(fac4)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*h4( ind + round( -1 ) ) - a3*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*3*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac4)*h4(ind);
        h4d(ind) = 4*h4d(ind);
        
    end
    for k4=k1+0+2:N
        ind = ind+1;
        % (k1,k1,k1,k4)
        h4d(ind) = a4*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*3*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac4)*h4( ind + round( -1 ) ) - a3*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*3*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac4)*h4(ind);
        h4d(ind) = 4*h4d(ind);
        
    end
    
    if (k1 <= N-1 && k1 <= N-1)
        ind = ind+1;
        k2 = k1;
        k3 = k1+1;
        % (k1,k1,k1+1,k1+1)
        h4d(ind) = a4*(fac1)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*2*(fac4)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*h4( ind + round( k3 - N - 2 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*4*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*2*(fac4)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac1)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac1)*h4(ind);
        h4d(ind) = 6*h4d(ind);
        
    end
    
    if (k1 <= N-1 && k1 <= N-2)
        % (k1,k1,k1+1,k1+2)
        ind = ind+1;
        h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac4)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac1)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac1)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac4)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    for k4=k1+1+2:N
        
        if (k1 <= N-1)
            % (k1,k1,k1+1,k4)
            ind = ind+1;
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac4)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac4)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
    end
    for k3=k1+0+2:N
        ind = ind+1;
        % (k1,k1,k3,k3)
        h4d(ind) = a4*(fac1)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac1)*h4( ind + round( k3 - N - 2 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*4*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*2*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac1)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac1)*h4(ind);
        h4d(ind) = 6*h4d(ind);
        
        
        if (k3 <= N-1)
            % (k1,k1,k3,k3+1)
            ind = ind+1;
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac1)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac1)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
        for k4=k3+0+2:N
            ind = ind+1;
            % (k1,k1,k3,k4)
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
    end
    
    if (k1 <= N-1 && k1 <= N-1 && k1 <= N-1)
        k2 = k1+1;
        k3 = k1+1;        
        % (k1,k1+1,k1+1,k1+1)
        ind = ind+1;
        h4d(ind) = a4*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*3*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*3*(fac1)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac4)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac4)*h4(ind);
        h4d(ind) = 4*h4d(ind);
        
    end
    
    if (k1 <= N-1 && k1 <= N-1 && k1 <= N-2)
        % (k1,k1+1,k1+1,k1+2)
        ind = ind+1;
        h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac1)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac4)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac4)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    for k4=k1+1+2:N
        
        if (k1 <= N-1 && k1 <= N-1)
            % (k1,k1+1,k1+1,k4)
            ind = ind+1;
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
    end
    
    if (k1 <= N-1 && k1 <= N-2 && k1 <= N-2)
        k2 = k1+1;
        k3 = k1+2;            
        % (k1,k1+1,k1+2,k1+2)
        ind = ind+1;
        h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac1)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac4)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac4)*h4( ind + round( k3 - N - 2 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*2*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac1)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
        h4d(ind) = 12*h4d(ind);
        
    end
    
    if (k1 <= N-1 && k1 <= N-2 && k1 <= N-3)
        % (k1,k1+1,k1+2,k1+3)
        ind = ind+1;
        h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac1)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
        h4d(ind) = 24*h4d(ind);
        
    end
    for k4=k1+2+2:N
        
        if (k1 <= N-1 && k1 <= N-2)
            % (k1,k1+1,k1+2,k4)
            ind = ind+1;
            h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
            h4d(ind) = 24*h4d(ind);
            
        end
    end
    for k3=k1+1+2:N
        
        if (k1 <= N-1)
            % (k1,k1+1,k3,k3)
            ind = ind+1;
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac1)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*2*(fac3)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac1)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
        
        if (k1 <= N-1 && k3 <= N-1)
            % (k1,k1+1,k3,k3+1)
            ind = ind+1;
            h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac1)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
            h4d(ind) = 24*h4d(ind);
            
        end
        for k4=k3+0+2:N
            
            if (k1 <= N-1)
                % (k1,k1+1,k3,k4)
                ind = ind+1;
                h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
                h4d(ind) = 24*h4d(ind);
                
            end
        end
    end
    for k2=k1+0+2:N
        ind = ind+1;
        % (k1,k2,k2,k2)
        k3 = k2;
        h4d(ind) = a4*(fac4)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac4)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*3*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac4)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac4)*h4(ind);
        h4d(ind) = 4*h4d(ind);
        
        
        if (k2 <= N-1)
            % (k1,k2,k2,k2+1)
            ind = ind+1;
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac4)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac4)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
        for k4=k2+0+2:N
            ind = ind+1;
            % (k1,k2,k2,k4)
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*2*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
        
        if (k2 <= N-1 && k2 <= N-1)
            % (k1,k2,k2+1,k2+1)
            k3 = k2+1;
            ind = ind+1;
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac4)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac4)*h4( ind + round( k3 - N - 2 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*2*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
        end
        
        if (k2 <= N-1 && k2 <= N-2)
            % (k1,k2,k2+1,k2+2)
            ind = ind+1;
            h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
            h4d(ind) = 24*h4d(ind);
            
        end
        for k4=k2+1+2:N
            
            if (k2 <= N-1)
                % (k1,k2,k2+1,k4)
                ind = ind+1;
                h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac2)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
                h4d(ind) = 24*h4d(ind);
                
            end
        end
        for k3=k2+0+2:N
            ind = ind+1;
            % (k1,k2,k3,k3)
            h4d(ind) = a4*(fac2)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac2)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac2)*h4( ind + round( k3 - N - 2 ) ) - a3*2*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*2*(fac3)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac2)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac2)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac2)*h4(ind);
            h4d(ind) = 12*h4d(ind);
            
            
            if (k3 <= N-1)
                % (k1,k2,k3,k3+1)
                ind = ind+1;
                h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac2)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac2)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac2)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac2)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
                h4d(ind) = 24*h4d(ind);
                
            end
            for k4=k3+0+2:N
                ind = ind+1;
                % (k1,k2,k3,k4)
                h4d(ind) = a4*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 4 ) ) - a3*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 3 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 3 ) ) + a2*(fac3)*h4( ind + round( k3 - N - 2 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - (c3N)/2 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( -((N - k1 + 4)*(c2N - 2*k1 - c2N*k1 + NA2 + k1^2 + 3))/6 ) ) - a1*(fac3)*h4( ind + round( -1 ) ) - a3*(fac3)*h4( ind + round( (11*k1)/6 - c13Nd3 + (3*k2)/2 + k3 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 3 ) ) + a2*(fac3)*h4( ind + round( (3*k2)/2 - c5Nd2 + k3 + N*k2 - NA2/2 - k2^2/2 - 2 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - c17Nd6 + k3 + c2N*k1 - (N*k1^2)/2 + (NA2*k1)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - 2 ) ) - a1*(fac3)*h4( ind + round( k3 - N - 1 ) ) + a2*(fac3)*h4( ind + round( (11*k1)/6 - (c10N)/3 + (3*k2)/2 + c2N*k1 + N*k2 - (N*k1^2)/2 - NA2 - NA3/6 - k1^2 + k1^3/6 - k2^2/2 + (NA2*(k1 - 1))/2 - 2 ) ) - a1*(fac3)*h4( ind + round( -((N - k2 + 1)*(N - k2 + 2))/2 ) ) - a1*(fac3)*h4( ind + round( -((N - k1 + 1)*(N - k1 + 2)*(N - k1 + 3))/6 ) ) + (fac3)*h4(ind);
                h4d(ind) = 24*h4d(ind);
                
            end
        end
    end
end

