% Integrals_PREC_N.m is a routine for constructing integrals which are used
% for evaluation of conditional expectations for the multi-dimensional
% case; it computes integrals of the form I=E(exp(l*e)), l being a N-dimensional 
% vector, for ordinary polynomials of degrees from one to five and a normally 
% distributed multivariate random vector of shocks with mean 0 and 
% variance-covariance matrix vcv; see "How to Solve Dynamic Stochastic Models 
% Computing Expectations Just Once" by Kenneth L. Judd, Lilia Maliar, Serguei 
% Maliar and Inna Tsener, (2016), Quantitative Economics, forthcoming, 
% (henceforth, JMMT, 2016).
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Inputs:  "n1"           is a number of endogenous state variables 
%          "n2"           is a number of exogenous state variables
%          "vcv"          is a variance-covariance matrix of shocks; 
%                         n2-by-n2
%          "D"            is a degree of the polynomial for which the sets
%                         of precomputed integrals must be constructed; (can
%                         be 1,2,3,4 or 5)
%
% Output: "integrals_fs"  is a vector of precomputed integrals for a complete 
%                         polynomial of the given degree 
% -------------------------------------------------------------------------
% Copyright © 2016 by Lilia Maliar, Serguei Maliar and Inna Tsener. All 
% rights reserved. The code may be used, modified and redistributed under 
% the terms provided in the file "License_Agreement.txt".
% -------------------------------------------------------------------------

function integrals_fs = Integrals_PREC_N(n1,n2,vcv,D)

% Note that "Integrals_PREC_N.m" is a function that must be used together 
% with "Ord_Polynomial_N.m". The input of this polynomial routine must 
% be in the form s = [k, z], where k and z are n1- and n2-dimensional vectors, 
% respectively. This ordering of state variables ensures that precomputed 
% integrals match the corresponding polynomial terms in the approximating 
% function of degree D

    n = n1+n2;             % Compute the total number of state variables
    
    % 1. The vector of precomputed integrals for first-degree polynomial 
    % (the default option)
    % -------------------------------------------------------------------
    
    for i = 1:n2            % Construct integrals for the coefficients that 
                            % multiply powers of z
        l = zeros(1,n2);    % Initial values of powers of z; these are zeros
        l(1,i) = 1;         % For the first-degree polynomial each country's 
                            % specific shock is entering the polynomial
                            % with power 1
        zexp(1,i) = exp(l*vcv*l'/2);
                            % Compute integrals exactly; see equaiton (53) 
                            % in JMMT (2016)
    end
    
    integrals_fs = [1 ones(1,n1) zexp];  
                        % The vector of precomputed integrals includes
                        % n1+1 ones (i.e., a constant and the coefficients  
                        % that multiply powers of k) and n2 integrals for 
                        % the coefficients that multiply powers of z
                                  
                                    
    i = n+1;     % Index number for precomputed integrals; precomputed integrals 
                 % for first-degree polynomial are indexed from 1 to n+1, 
                 % and subsequent integrals will be indexed from n+2 and on
    
    % 2. The vector of precomputed integrals for second-degree polynomial 
    % ------------------------------------------------------------------- 
    if D == 2 

       for j1 = 1:n             % For each state variable...
           for j2 = j1:n        % take another state variable and compute 
                                % exactly the integral that multiplies the
                                % coefficient on the product of these two 
                                % variables in the polynomial
               
                i = i+1;        % Increase index number by one
                
                if (j2>n1)&&(j1<=n1)                  
                    % Compute integrals for the coefficients that multiply 
                    % ONLY products of endogenous and exogenous state 
                    % variables, for example: k1*z1, k2*z2,...                   
                    l = zeros(1,n2);    % Initial values of powers of z; 
                                        % these are zeros
                    l(1,j2-n1) = 1;     % Each country's specific shock is 
                                        % entering the polynomial with power 1
                    zexp = exp(l*vcv*l'/2);
                                        % Compute integrals exactly 
                    integrals_fs(:,i) = zexp;
                                        % Store precomputed integrals

                elseif (j2>n1)&&(j1>n1)
                    % Compute integrals for the coefficients that multiply 
                    % ONLY products of exogenous state variables; for
                    % example, z1*z1, z1*z2, z2*z2,...
                    if j1==j2
                                        % Compute the squares of exogenous 
                                        % state variables
                        l = zeros(1,n2);
                        l(1,j2-n1) = 2; % Each country's specific shock is 
                                        % entering polynomial with power 2
                        zexp = exp(l*vcv*l'/2); 
                        integrals_fs(:,i) = zexp;
                    else    
                                        % Compute interactions of the exogenous 
                                        % state variables
                        l = zeros(1,n2);
                        l(1,j2-n1) = 1; % Each country's specific shock is
                        l(1,j1-n1) = 1;  % entering polynomial with power 1
                        zexp = exp(l*vcv*l'/2);
                        integrals_fs(:,i) = zexp;
                    end
                else
                    % Compute integrals for the coefficients that multiply 
                    % ONLY products of endogenous state variables; for
                    % example, k1*k1, k2*k2, k1*k2,...
                    integrals_fs(:,i) = 1; 
                end
           end
       end
 
    % 3.  The vector of precomputed integrals for third-degree polynomial 
    % -------------------------------------------------------------------    
    elseif D == 3                
        for j1 = 1:n
            for j2 = j1:n
                i = i+1;
                if (j2>n1)&&(j1<=n1)
                    
                l = zeros(1,n2);
                l(1,j2-n1) = 1; 
                zexp = exp(l*vcv*l'/2);
                integrals_fs(:,i) = zexp;

                elseif (j2>n1)&&(j1>n1)

                    if j1==j2
                    l = zeros(1,n2);
                    l(1,j2-n1) = 2;
                    zexp = exp(l*vcv*l'/2);
                    integrals_fs(:,i) = zexp;
                    else    
                    l = zeros(1,n2);
                    l(1,j2-n1) = 1;
                    l(1,j1-n1) = 1;
                    zexp = exp(l*vcv*l'/2);
                    integrals_fs(:,i) = zexp;
                    end
                else
                    integrals_fs(:,i) = 1;
                end
                
                for j3 = j2:n
                    i = i+1;
                    if (j1<=n1)&&(j2<=n1)
                       if (j3<=n1)
                       integrals_fs(:,i) = 1;
                       else
                       l = zeros(1,n2);
                       l(1,j3-n1) = 1; 
                       zexp = exp(l*vcv*l'/2); 
                       integrals_fs(:,i) = zexp;
                       end
                    elseif  (j1<=n1)&&(j2>n1)  
                       if (j3==j2)
                       l = zeros(1,n2);
                       l(1,j2-n1) = 2;
                       zexp = exp(l*vcv*l'/2);
                       integrals_fs(:,i) = zexp;
                       else
                       l = zeros(1,n2);
                       l(1,j2-n1) = 1;
                       l(1,j3-n1) = 1;
                       zexp = exp(l*vcv*l'/2);
                       integrals_fs(:,i) = zexp;  
                       end
                    elseif  (j1>n1)&&(j2>n1) 
                        if (j1==j2)
                           if (j2==j3)
                           l = zeros(1,n2);
                           l(1,j3-n1) = 3;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j2-n1) = 2;
                           l(1,j3-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           end
                        else
                           if (j2==j3)
                           l = zeros(1,n2);
                           l(1,j1-n1) = 1;
                           l(1,j3-n1) = 2;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j1-n1) = 1;
                           l(1,j2-n1) = 1;
                           l(1,j3-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           end
                       end
                    end
                end
            end
        end

    

    % 4.  The vector of precomputed integrals for fourth-degree polynomial 
    % ------------------------------------------------------------------------    
    elseif D == 4
        for j1 = 1:n
            for j2 = j1:n
                i = i+1;
                
                if (j2>n1)&&(j1<=n1) 

                l = zeros(1,n2);
                l(1,j2-n1) = 1; 
                zexp = exp(l*vcv*l'/2);
                integrals_fs(:,i) = zexp;

                elseif (j2>n1)&&(j1>n1)

                    if j1==j2
                    l = zeros(1,n2);
                    l(1,j2-n1) = 2;
                    zexp = exp(l*vcv*l'/2);
                    integrals_fs(:,i) = zexp;
                    else    
                    l = zeros(1,n2);
                    l(1,j2-n1) = 1;
                    l(1,j1-n1) = 1;
                    zexp = exp(l*vcv*l'/2);
                    integrals_fs(:,i) = zexp;
                    end
                else
                    integrals_fs(:,i) = 1;
                end
                
                for j3 = j2:n
                    i = i+1;
                    
                   if (j1<=n1)&&(j2<=n1)
                       if (j3<=n1)
                       integrals_fs(:,i) = 1;
                       else
                       l = zeros(1,n2);
                       l(1,j3-n1) = 1; 
                       zexp = exp(l*vcv*l'/2); 
                       integrals_fs(:,i) = zexp;
                       end
                   elseif  (j1<=n1)&&(j2>n1)
                       if (j3==j2)
                       l = zeros(1,n2);
                       l(1,j2-n1) = 2;
                       zexp = exp(l*vcv*l'/2);
                       integrals_fs(:,i) = zexp;
                       else
                       l = zeros(1,n2);
                       l(1,j2-n1) = 1;
                       l(1,j3-n1) = 1;
                       zexp = exp(l*vcv*l'/2);
                       integrals_fs(:,i) = zexp;  
                       end
                   elseif  (j1>n1)&&(j2>n1)   
                       
                       if (j1==j2)
                           if (j2==j3)
                           l = zeros(1,n2);
                           l(1,j3-n1) = 3;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j2-n1) = 2;
                           l(1,j3-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           end
                       else
                           if (j2==j3)
                           l = zeros(1,n2);
                           l(1,j1-n1) = 1;
                           l(1,j3-n1) = 2;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j1-n1) = 1;
                           l(1,j2-n1) = 1;
                           l(1,j3-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           end
                       end
                       
                   end
                       for j4 = j3:n
                       i = i+1;
                       if (j1<=n1)&&(j2<=n1)&&(j3<=n1)
                           if (j4<=n1)
                           integrals_fs(:,i) = 1;
                           else
                           l = zeros(1,n2);
                           l(1,j4-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;    
                           end
                       elseif (j1<=n1)&&(j2<=n1)&&(j3>n1)
                           if (j3==j4)
                           l = zeros(1,n2);
                           l(1,j4-n1) = 2;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j3-n1) = 1;
                           l(1,j4-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp; 
                           end
                       elseif (j1<=n1)&&(j2>n1)&&(j3>n1)
                            
                            if (j2==j3)
                                if j3==j4
                                l = zeros(1,n2);
                                l(1,j2-n1) = 3;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                else
                                l = zeros(1,n2);
                                l(1,j3-n1) = 2;
                                l(1,j4-n1) = 1;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                end
                                
                            else
                                if j3==j4
                                l = zeros(1,n2);
                                l(1,j2-n1) = 1;
                                l(1,j3-n1) = 2;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                else
                                l = zeros(1,n2);
                                l(1,j2-n1) = 1;
                                l(1,j3-n1) = 1;
                                l(1,j4-n1) = 1;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                end
                            end
                       elseif (j1>n1)&&(j2>n1)&&(j3>n1)
                            if j1==j2
                                if j2==j3
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j4-n1) = 4;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j3-n1) = 3;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                else
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j2-n1) = 2;
                                    l(1,j4-n1) = 2;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 2;
                                    l(1,j3-n1) = 1;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                end 
                            else
                                if j2==j3
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j4-n1) = 3;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j3-n1) = 2;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                else
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j2-n1) = 1;
                                    l(1,j4-n1) = 2;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j2-n1) = 1;
                                    l(1,j3-n1) = 1;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                end
                            end
                       end
                       end   
               end
            end
        end 
            
            
  
    % 5.  The vector of precomputed integrals for fifth-degree polynomial 
    % ------------------------------------------------------------------- 

    elseif D == 5         
       
        for j1 = 1:n
            for j2 = j1:n
                i = i+1;
                if (j2>n1)&&(j1<=n1) 

                l = zeros(1,n2);
                l(1,j2-n1) = 1; 
                zexp = exp(l*vcv*l'/2);
                integrals_fs(:,i) = zexp;

                elseif (j2>n1)&&(j1>n1)

                    if j1==j2
                    l = zeros(1,n2);
                    l(1,j2-n1) = 2;
                    zexp = exp(l*vcv*l'/2);
                    integrals_fs(:,i) = zexp;
                    else    
                    l = zeros(1,n2);
                    l(1,j2-n1) = 1;
                    l(1,j1-n1) = 1;
                    zexp = exp(l*vcv*l'/2);
                    integrals_fs(:,i) = zexp;
                    end
                else
                    integrals_fs(:,i) = 1;
                end
                
                for j3 = j2:n
                    i = i+1;

                    
                   if (j1<=n1)&&(j2<=n1)
                       if (j3<=n1)
                       integrals_fs(:,i) = 1;
                       else
                       l = zeros(1,n2);
                       l(1,j3-n1) = 1; 
                       zexp = exp(l*vcv*l'/2); 
                       integrals_fs(:,i) = zexp;
                       end
                   elseif  (j1<=n1)&&(j2>n1)  
                       if (j3==j2)
                       l = zeros(1,n2);
                       l(1,j2-n1) = 2;
                       zexp = exp(l*vcv*l'/2);
                       integrals_fs(:,i) = zexp;
                       else
                       l = zeros(1,n2);
                       l(1,j2-n1) = 1;
                       l(1,j3-n1) = 1;
                       zexp = exp(l*vcv*l'/2);
                       integrals_fs(:,i) = zexp;  
                       end
                   elseif  (j1>n1)&&(j2>n1)   
                       
                       if (j1==j2)
                           if (j2==j3)
                           l = zeros(1,n2);
                           l(1,j3-n1) = 3;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j2-n1) = 2;
                           l(1,j3-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           end
                       else
                           if (j2==j3)
                           l = zeros(1,n2);
                           l(1,j1-n1) = 1;
                           l(1,j3-n1) = 2;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j1-n1) = 1;
                           l(1,j2-n1) = 1;
                           l(1,j3-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           end
                       end
                       
                   end
                   
                    for j4 = j3:n
                        i = i+1;
                        if (j1<=n1)&&(j2<=n1)&&(j3<=n1)
                           if (j4<=n1)
                           integrals_fs(:,i) = 1;
                           else
                           l = zeros(1,n2);
                           l(1,j4-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;    
                           end
                        elseif (j1<=n1)&&(j2<=n1)&&(j3>n1)
                           if (j3==j4)
                           l = zeros(1,n2);
                           l(1,j4-n1) = 2;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp;
                           else
                           l = zeros(1,n2);
                           l(1,j3-n1) = 1;
                           l(1,j4-n1) = 1;
                           zexp = exp(l*vcv*l'/2);
                           integrals_fs(:,i) = zexp; 
                           end
                        elseif (j1<=n1)&&(j2>n1)&&(j3>n1)
                           if (j2==j3)
                                if j3==j4
                                l = zeros(1,n2);
                                l(1,j2-n1) = 3;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                else
                                l = zeros(1,n2);
                                l(1,j3-n1) = 2;
                                l(1,j4-n1) = 1;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                end
                           else
                                if j3==j4
                                l = zeros(1,n2);
                                l(1,j2-n1) = 1;
                                l(1,j3-n1) = 2;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                else
                                l = zeros(1,n2);
                                l(1,j2-n1) = 1;
                                l(1,j3-n1) = 1;
                                l(1,j4-n1) = 1;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                end
                           end
                        elseif (j1>n1)&&(j2>n1)&&(j3>n1)
                            if j1==j2
                                if j2==j3
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j4-n1) = 4;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j3-n1) = 3;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                else
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j2-n1) = 2;
                                    l(1,j4-n1) = 2;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 2;
                                    l(1,j3-n1) = 1;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                end 
                            else
                                if j2==j3
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j4-n1) = 3;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j3-n1) = 2;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                else
                                    if j3==j4
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j2-n1) = 1;
                                    l(1,j4-n1) = 2;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j1-n1) = 1;
                                    l(1,j2-n1) = 1;
                                    l(1,j3-n1) = 1;
                                    l(1,j4-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end
                                end
                            end
                       end
                        for j5 = j4:n
                            i = i+1;
                            if (j1<=n1)&&(j2<=n1)&&(j3<=n1)&&(j4<=n1)
                                if (j5<=n1)
                                integrals_fs(:,i) = 1;
                                else
                                l = zeros(1,n2);
                                l(1,j5-n1) = 1;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                end
                            elseif (j1<=n1)&&(j2<=n1)&&(j3<=n1)&&(j4>n1)
                                if (j4==j5)
                                l = zeros(1,n2);
                                l(1,j5-n1) = 2;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                else
                                l = zeros(1,n2);
                                l(1,j4-n1) = 1;
                                l(1,j5-n1) = 1;
                                zexp = exp(l*vcv*l'/2);
                                integrals_fs(:,i) = zexp; 
                                end
                            elseif (j1<=n1)&&(j2<=n1)&&(j3>n1)&&(j4>n1)
                                if j3==j4
                                    if j4==j5
                                    l = zeros(1,n2);
                                    l(1,j4-n1) = 3;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j4-n1) = 2;
                                    l(1,j5-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end   
                                else
                                    if j4==j5
                                    l = zeros(1,n2);
                                    l(1,j3-n1) = 1;
                                    l(1,j4-n1) = 2;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    else
                                    l = zeros(1,n2);
                                    l(1,j3-n1) = 1;
                                    l(1,j4-n1) = 1;
                                    l(1,j5-n1) = 1;
                                    zexp = exp(l*vcv*l'/2);
                                    integrals_fs(:,i) = zexp; 
                                    end   
                                end
                            elseif (j1<=n1)&&(j2>n1)&&(j3>n1)&&(j4>n1)    
                                if j2==j3
                                    if j3==j4
                                        if j4==j5
                                        l = zeros(1,n2);
                                        l(1,j3-n1) = 4;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp; 
                                        else
                                        l = zeros(1,n2);
                                        l(1,j3-n1) = 3;
                                        l(1,j5-n1) = 1;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp; 
                                        end
                                    else
                                        if j4==j5
                                        l = zeros(1,n2);
                                        l(1,j3-n1) = 2;
                                        l(1,j5-n1) = 2;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp; 
                                        else
                                        l = zeros(1,n2);
                                        l(1,j3-n1) = 2;
                                        l(1,j4-n1) = 1;
                                        l(1,j5-n1) = 1;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp;
                                        end
                                    end
                                else
                                    if j3==j4
                                        if j4==j5
                                        l = zeros(1,n2);
                                        l(1,j2-n1) = 1;
                                        l(1,j5-n1) = 3;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp;
                                        else
                                        l = zeros(1,n2);
                                        l(1,j2-n1) = 1;
                                        l(1,j4-n1) = 2;
                                        l(1,j5-n1) = 1;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp;
                                        end
                                    else
                                        if j4==j5
                                        l = zeros(1,n2);
                                        l(1,j2-n1) = 1;
                                        l(1,j3-n1) = 1;
                                        l(1,j5-n1) = 2;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp;
                                        else
                                        l = zeros(1,n2);
                                        l(1,j2-n1) = 1;
                                        l(1,j3-n1) = 1;
                                        l(1,j4-n1) = 1;
                                        l(1,j5-n1) = 1;
                                        zexp = exp(l*vcv*l'/2);
                                        integrals_fs(:,i) = zexp;
                                        end
                                    end
                                end
                            elseif (j1>n1)&&(j2>n1)&&(j3>n1)&&(j4>n1) 
                                if j1==j2
                                    if j2==j3
                                        if j3==j4
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 5;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 4;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        else
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 3;
                                            l(1,j4-n1) = 2;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 3;
                                            l(1,j4-n1) = 1;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        end
                                    else
                                        if j3==j4
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 2;
                                            l(1,j4-n1) = 3;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 2;
                                            l(1,j4-n1) = 2;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        else
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 2;
                                            l(1,j3-n1) = 1;
                                            l(1,j4-n1) = 2;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 2;
                                            l(1,j3-n1) = 1;
                                            l(1,j4-n1) = 1;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        end
                                    end
                                else
                                    if j2==j3
                                        if j3==j4
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j2-n1) = 4;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j3-n1) = 3;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        else
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j2-n1) = 2;
                                            l(1,j5-n1) = 2;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j3-n1) = 2;
                                            l(1,j4-n1) = 1;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        end
                                   else
                                        if j3==j4
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j2-n1) = 1;
                                            l(1,j5-n1) = 3;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j2-n1) = 1;
                                            l(1,j4-n1) = 2;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        else
                                            if j4==j5
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j2-n1) = 1;
                                            l(1,j3-n1) = 1;
                                            l(1,j5-n1) = 2;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            else
                                            l = zeros(1,n2);
                                            l(1,j1-n1) = 1;
                                            l(1,j2-n1) = 1;
                                            l(1,j3-n1) = 1;
                                            l(1,j4-n1) = 1;
                                            l(1,j5-n1) = 1;
                                            zexp = exp(l*vcv*l'/2);
                                            integrals_fs(:,i) = zexp;
                                            end
                                        end
                                    end
                                end
                            end
                       end
                    end   
                end
            end
        end    
   end                         % end of "if"/"elseif"   