function[p]=gau(a,b,sigma,x0,nn)
t=linspace(a,b,nn);
m=1/(sqrt(2*pi)*sigma);
for i=1:nn,
   ff(i)=m*exp(-(t(i)-x0)^2/(2*sigma*sigma));
end
h=(b-a)/(nn-1);
u=ff(2:2:nn-1);
vv=ff(3:2:nn-2);ss=sum(u);
v=sum(vv);
p=(h/3)*(ff(1)+ff(nn)+4*ss+2*v);