function [Ke,Fe]=symbolicKF(num,D)
h=1.0/num;
syms x y z
x0=0; y0=0; z0=0;
N(1) = (h+x0-x) * (h+y0-y) * (h+z0-z) / (h^3);
N(2) = (x-x0) * (h+y0-y) * (h+z0-z) / (h^3);
N(3) = (h+x0-x) * (y-y0) * (h+z0-z) / (h^3);
N(4) = (x-x0) * (y-y0) * (h+z0-z) / (h^3);
N(5) = (h+x0-x) * (h+y0-y) * (z-z0) / (h^3);
N(6) = (x-x0) * (h+y0-y) * (z-z0) / (h^3);
N(7) = (h+x0-x) * (y-y0) * (z-z0) / (h^3);
N(8) = (x-x0) * (y-y0) * (z-z0) / (h^3);

for j=1:8
    for jj=1:3
        Dif(3*jj-2,j*3+jj-3) = diff(N(j),x);
        Dif(3*jj-1,j*3+jj-3) = diff(N(j),y);
        Dif(3*jj,  j*3+jj-3) = diff(N(j),z);
    end
end
B9x24 = Dif;
for j=1:24
    B(1,j) = B9x24(1,j);
    B(2,j) = B9x24(5,j);
    B(3,j) = B9x24(9,j);
    B(4,j) = B9x24(2,j) + B9x24(4,j);
    B(5,j) = B9x24(6,j) + B9x24(8,j);
    B(6,j) = B9x24(3,j) + B9x24(7,j);
end

K=B.'*D*B;
Ke=zeros(24,24);
for i=1:24
    for j=i:24
        Ke(i,j)=int(int(int(K(i,j),'x',0,h),'y',0,h),'z',0,h);
    end
end
Ke=Ke+Ke.'-diag(diag(Ke));

F=B.'*D;
Fe=zeros(24,6);
for i=1:24
    for j=1:6
        Fe(i,j)=int(int(int(F(i,j),'x',0,h),'y',0,h),'z',0,h);
    end
end
end