function t=simpson(a,b,n)
  h=(b-a)/n;
  x=a:h:b;
  y=f(x);
  t=y(1)+y(n+1);
  for i=2:2:n
    t=t+4*y(i);
  endfor
  for j=3:2:n-1
    t=t+2*y(j);
  endfor
    t=t*h/3;
endfunction