function sierpinski_3D
    function tetraeder(x,y,z,inv,size)
        points=[1,0,-1/sqrt(2);
                0,1,1/sqrt(2);
                -1,0,-1/sqrt(2);
                0,-1,1/sqrt(2);
                0,1,1/sqrt(2);
                1,0,-1/sqrt(2);
                0,-1,1/sqrt(2);
                -1,0,-1/sqrt(2);
                1,0,-1/sqrt(2);]/2;
        if(inv==1)
            points(:,3)=-1*points(:,3);
        end
        plot3(points(:,1).*size + x, points(:,2).*size + y, points(:,3).*size + z);
    end
    figure(1);
    hold on;
    tetraeder(0,0,0,0,1);
    a=[1/2,1/2,0];
    b=[-1/2,-1/2,0];
    %tetraeder(m2/2,m2/2,m2/2,0,1/2);
    function iterate(x,y,z,size,itnum)
        m2=sqrt(sum((b-a).*(b-a)));
        size=size/2;
        itnum=itnum+1;
        tetraeder(-0.25*size+x,0+y,size*-1*m2/8+z,0,size*1/2);
        tetraeder(size*0.25+x,0+y,size*-1*m2/8+z,0,size*1/2);
        tetraeder(0+x,size*-1*0.25+y,size*m2/8+z,0,size*1/2);
        tetraeder(0+x,size*0.25+y,size*m2/8+z,0,size*1/2);
        if(itnum<4)
            iterate(x+-0.25*size,y,z-1*size*m2/8,size,itnum);
            iterate(x+0.25*size,y,z-1*size*m2/8,size,itnum);
            iterate(x,y-0.25*size,z+m2*size/8,size,itnum);
            iterate(x,y+0.25*size,z+m2*size/8,size,itnum);
        end
    end
    m2=sqrt(sum((b-a).*(b-a)));
    iterate(-0.25,0,-m2/8,1,1);
    iterate(0.25,0,-m2/8,1,1);
    iterate(0,-0.25,m2/8,1,1);
    iterate(0,0.25,m2/8,1,1);
    grid on;
    xlabel('x');
    ylabel('y');
    zlabel('z');
end