how to triple Integrate this?

clear
clc
D = 200000;

syms xi eta zeta

N(1) = (1/8)*(1-xi)*(1-eta)*(1-zeta);
N(2) = (1/8)*(1+xi)*(1-eta)*(1-zeta);
N(3) = (1/8)*(1+xi)*(1+eta)*(1-zeta);
N(4) = (1/8)*(1-xi)*(1+eta)*(1-zeta);
N(5) = (1/8)*(1-xi)*(1-eta)*(1+zeta);
N(6) = (1/8)*(1+xi)*(1-eta)*(1+zeta);
N(7) = (1/8)*(1+xi)*(1+eta)*(1+zeta);
N(8) = (1/8)*(1-xi)*(1+eta)*(1+zeta);

for i = 1:8
    N_xi(i) = diff(N(i),xi);
    N_eta(i) = diff(N(i),eta);
    N_zeta(i) = diff(N(i),zeta);
end

% Nodes   xi  eta zeta
%  1      -1  -1  -1
%  2       1  -1  -1
%  3       1   1  -1
%  4      -1   1  -1
%  5      -1  -1   1
%  6       1  -1   1
%  7       1   1   1
%  8      -1   1   1


% x1_xi   = (N_xi(1)*-1) + (N_xi(2)*1) + (N_xi(3)*1) + (N_xi(4)*-1) + (N_xi(5)*-1) + (N_xi(6)*1) + (N_xi(7)*1) + (N_xi(8)*-1);
% x1_eta  = (N_eta(1)*-1) + (N_eta(2)*1) + (N_eta(3)*1) + (N_eta(4)*-1) + (N_eta(5)*-1) + (N_eta(6)*1) + (N_eta(7)*1) + (N_eta(8)*-1);
% x1_zeta = (N_zeta(1)*-1) + (N_zeta(2)*1) + (N_zeta(3)*1) + (N_zeta(4)*-1) + (N_zeta(5)*-1) + (N_zeta(6)*1) + (N_zeta(7)*1) + (N_zeta(8)*-1);
% x2_xi   = (N_xi(1)*-1) + (N_xi(2)*-1) + (N_xi(3)*1) + (N_xi(4)*1) + (N_xi(5)*-1) + (N_xi(6)*-1) + (N_xi(7)*1) + (N_xi(8)*1);
% x2_eta  = (N_eta(1)*-1) + (N_eta(2)*-1) + (N_eta(3)*1) + (N_eta(4)*1) + (N_eta(5)*-1) + (N_eta(6)*-1) + (N_eta(7)*1) + (N_eta(8)*1);
% x2_zeta = (N_zeta(1)*-1) + (N_zeta(2)*-1) + (N_zeta(3)*1) + (N_zeta(4)*1) + (N_zeta(5)*-1) + (N_zeta(6)*-1) + (N_zeta(7)*1) + (N_zeta(8)*1);
% x3_xi   = (N_xi(1)*-1) + (N_xi(2)*-1) + (N_xi(3)*-1) + (N_xi(4)*-1) + (N_xi(5)*1) + (N_xi(6)*1) + (N_xi(7)*1) + (N_xi(8)*1);
% x3_eta  = (N_eta(1)*-1) + (N_eta(2)*-1) + (N_eta(3)*-1) + (N_eta(4)*-1) + (N_eta(5)*1) + (N_eta(6)*1) + (N_eta(7)*1) + (N_eta(8)*1);
% x3_zeta = (N_zeta(1)*-1) + (N_zeta(2)*-1) + (N_zeta(3)*-1) + (N_zeta(4)*-1) + (N_zeta(5)*1) + (N_zeta(6)*1) + (N_zeta(7)*1) + (N_zeta(8)*1);
% 
% J = [x1_xi x1_eta x1_zeta;x2_xi x2_eta x2_zeta;x3_xi x3_eta x3_zeta];

J = [1 0 0;0 1 0;0 0 1];
Inv_J = inv(J);
Det_J = det(J);

N_x = [Inv_J(1,1) * N_xi; Inv_J(2,2) * N_eta; Inv_J(3,3) * N_zeta];

N_x = [N_x(1,1) 0 0 N_x(1,2) 0 0 N_x(1,3) 0 0 N_x(1,4) 0 0 N_x(1,5) 0 0 N_x(1,6) 0 0 N_x(1,7) 0 0 N_x(1,8) 0 0;...
       0 N_x(2,1) 0 0 N_x(2,2) 0 0 N_x(2,3) 0 0 N_x(2,4) 0 0 N_x(2,5) 0 0 N_x(2,6) 0 0 N_x(2,7) 0 0 N_x(2,8) 0;...
       0 0 N_x(3,1) 0 0 N_x(3,2) 0 0 N_x(3,3) 0 0 N_x(3,4) 0 0 N_x(3,5) 0 0 N_x(3,6) 0 0 N_x(3,7) 0 0 N_x(3,8)];

   
K_ma = N_x.' * N_x;
   


for i = 1:24
    for j = 1:24
        K_mat(i,j) = @(xi,eta,zeta) K_ma(i,j);
        K(i,j) = integral3(K_ma(i,j),-1,1,-1,1,-1,1);
    end
end

Answer 1

[Walter Roberson]

clear
clc
D = 200000;

syms xi eta zeta

N(1) = (1/8)*(1-xi)*(1-eta)*(1-zeta);
N(2) = (1/8)*(1+xi)*(1-eta)*(1-zeta);
N(3) = (1/8)*(1+xi)*(1+eta)*(1-zeta);
N(4) = (1/8)*(1-xi)*(1+eta)*(1-zeta);
N(5) = (1/8)*(1-xi)*(1-eta)*(1+zeta);
N(6) = (1/8)*(1+xi)*(1-eta)*(1+zeta);
N(7) = (1/8)*(1+xi)*(1+eta)*(1+zeta);
N(8) = (1/8)*(1-xi)*(1+eta)*(1+zeta);

for i = 1:8
    N_xi(i) = diff(N(i),xi);
    N_eta(i) = diff(N(i),eta);
    N_zeta(i) = diff(N(i),zeta);
end

% Nodes   xi  eta zeta
%  1      -1  -1  -1
%  2       1  -1  -1
%  3       1   1  -1
%  4      -1   1  -1
%  5      -1  -1   1
%  6       1  -1   1
%  7       1   1   1
%  8      -1   1   1


% x1_xi   = (N_xi(1)*-1) + (N_xi(2)*1) + (N_xi(3)*1) + (N_xi(4)*-1) + (N_xi(5)*-1) + (N_xi(6)*1) + (N_xi(7)*1) + (N_xi(8)*-1);
% x1_eta  = (N_eta(1)*-1) + (N_eta(2)*1) + (N_eta(3)*1) + (N_eta(4)*-1) + (N_eta(5)*-1) + (N_eta(6)*1) + (N_eta(7)*1) + (N_eta(8)*-1);
% x1_zeta = (N_zeta(1)*-1) + (N_zeta(2)*1) + (N_zeta(3)*1) + (N_zeta(4)*-1) + (N_zeta(5)*-1) + (N_zeta(6)*1) + (N_zeta(7)*1) + (N_zeta(8)*-1);
% x2_xi   = (N_xi(1)*-1) + (N_xi(2)*-1) + (N_xi(3)*1) + (N_xi(4)*1) + (N_xi(5)*-1) + (N_xi(6)*-1) + (N_xi(7)*1) + (N_xi(8)*1);
% x2_eta  = (N_eta(1)*-1) + (N_eta(2)*-1) + (N_eta(3)*1) + (N_eta(4)*1) + (N_eta(5)*-1) + (N_eta(6)*-1) + (N_eta(7)*1) + (N_eta(8)*1);
% x2_zeta = (N_zeta(1)*-1) + (N_zeta(2)*-1) + (N_zeta(3)*1) + (N_zeta(4)*1) + (N_zeta(5)*-1) + (N_zeta(6)*-1) + (N_zeta(7)*1) + (N_zeta(8)*1);
% x3_xi   = (N_xi(1)*-1) + (N_xi(2)*-1) + (N_xi(3)*-1) + (N_xi(4)*-1) + (N_xi(5)*1) + (N_xi(6)*1) + (N_xi(7)*1) + (N_xi(8)*1);
% x3_eta  = (N_eta(1)*-1) + (N_eta(2)*-1) + (N_eta(3)*-1) + (N_eta(4)*-1) + (N_eta(5)*1) + (N_eta(6)*1) + (N_eta(7)*1) + (N_eta(8)*1);
% x3_zeta = (N_zeta(1)*-1) + (N_zeta(2)*-1) + (N_zeta(3)*-1) + (N_zeta(4)*-1) + (N_zeta(5)*1) + (N_zeta(6)*1) + (N_zeta(7)*1) + (N_zeta(8)*1);
% 
% J = [x1_xi x1_eta x1_zeta;x2_xi x2_eta x2_zeta;x3_xi x3_eta x3_zeta];

J = [1 0 0;0 1 0;0 0 1];
Inv_J = inv(J);
Det_J = det(J);

N_x = [Inv_J(1,1) * N_xi; Inv_J(2,2) * N_eta; Inv_J(3,3) * N_zeta];

N_x = [N_x(1,1) 0 0 N_x(1,2) 0 0 N_x(1,3) 0 0 N_x(1,4) 0 0 N_x(1,5) 0 0 N_x(1,6) 0 0 N_x(1,7) 0 0 N_x(1,8) 0 0;...
       0 N_x(2,1) 0 0 N_x(2,2) 0 0 N_x(2,3) 0 0 N_x(2,4) 0 0 N_x(2,5) 0 0 N_x(2,6) 0 0 N_x(2,7) 0 0 N_x(2,8) 0;...
       0 0 N_x(3,1) 0 0 N_x(3,2) 0 0 N_x(3,3) 0 0 N_x(3,4) 0 0 N_x(3,5) 0 0 N_x(3,6) 0 0 N_x(3,7) 0 0 N_x(3,8)];

K_ma = N_x.' * N_x;
  
K = int(int(int(K_ma,xi,-1,1),eta,-1,1),zeta,-1,1)

K =