tetra_int.c running 
Tetrahedron P1(1,1,1)     P2(3,2.5,2) 
            P3(2,0.5,1.5) P4(0.5,1.5,3.5) 
f(x,y,z)=x^3+2y^2+x*y*z 
Volume by determinant = 1.14583 
bounds xmin=0.5, xmax=3, ymin=0.5, ymax=2.5, zmin=1, zmax=3.5 
bound volume = 12.5 
f1(x,y,z)=x+2*y+3*z+4  np=1 
computed integral=16.4714, exact integral=16.4714, error=6.39488e-14 
f1(x,y,z)=x+2*y+3*z+4  np=8 
computed integral=16.4714, exact integral=16.4714, error=6.75016e-14 
f1(x,y,z)=x+2*y+3*z+4  np=64 
computed integral=16.4714, exact integral=16.4714, error=6.03961e-14 
f1(x,y,z)=x+2*y+3*z+4  np=512 
computed integral=16.4714, exact integral=16.4714, error=7.4607e-14 

f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=1 
computed integral=126.597, exact integral=127.947, error=-1.34994 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=8 
computed integral=127.434, exact integral=127.947, error=-0.51294 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=64 
computed integral=127.806, exact integral=127.947, error=-0.140425 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=512 
computed integral=127.909, exact integral=127.947, error=-0.0378332 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=4096 
computed integral=127.936, exact integral=127.947, error=-0.0101954 

f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=1 
computed integral=691.289, exact integral=719.931, error=-28.6418 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=8 
computed integral=709.867, exact integral=719.931, error=-10.0639 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=64 
computed integral=717.195, exact integral=719.931, error=-2.73534 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=512 
computed integral=719.197, exact integral=719.931, error=-0.733985 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=4096 
computed integral=719.733, exact integral=719.931, error=-0.197288 

fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=1 
computed integral=170.057, exact integral=229.535, error=-59.4779 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=8 
computed integral=208.899, exact integral=229.535, error=-20.6358 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=64 
computed integral=223.886, exact integral=229.535, error=-5.64878 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=512 
computed integral=227.987, exact integral=229.535, error=-1.54801 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=4096 
computed integral=229.115, exact integral=229.535, error=-0.419497 

Using FEM Handbook method 
f1(x,y,z)=x+2y+3z+4 np=4 
computed integral=16.4714, exact integral=16.4714, error=6.39488e-14 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=4 
computed integral=127.947, exact integral=127.947, error=-3.89499e-07 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=4 
computed integral=719.943, exact integral=719.931, error=0.011892 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=4 
computed integral=227.777, exact integral=229.535, error=-1.75792 

Using tetraquad.m  method 
f1(x,y,z)=x+2y+3z+4 np=27 
computed integral=16.4714, exact integral=16.4714, error=-3.54028e-11 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=27 
computed integral=127.947, exact integral=127.947, error=-4.16971e-07 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=27 
computed integral=719.931, exact integral=719.931, error=-3.59394e-05 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=27 
computed integral=229.468, exact integral=229.535, error=-0.0670278