f=inline('exp(-3*x.^2)','x');
a=0; b=1; n=10; h=(b-a)/n; 
x=linspace(a+h/2,b-h/2,n);
w=ones(size(x));
y=f(x);
M=h*dot(y,w);

x=linspace(a,b,n+1);
y=f(x);
w=[1,2*ones(1,n-1),1];
T=h/2*dot(y,w);


w=[1,(4-(1+(-1).^(1:n-1))),1]
S=h/3*dot(y,w);

printf('M=%0.5f, T=%0.5f\n',M,T);
printf('S=%0.5f, (2*M+T)/3=%0.5f\n',S,(2*M+T)/3);
 
