function y=gauss4pts(fh,a,b)
%assign sample points and weights
y4=[sqrt((3-4*sqrt(0.3))/7),sqrt((3+4*sqrt(0.3))/7)];
w4=[0.5+sqrt(10/3)/12,0.5-sqrt(10/3)/12];
c1=(b-a)/2;
c2=(a+b)/2;
s4=0;
for i=1:2
    s4=s4+w4(i)*(fh(c2+c1*y4(i))+fh(c2-c1*y4(i)));
end
y=c1*s4;