syms x y z x0 y0 z0 h
% 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