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:n
    if rem(i,2)=0
      t=t+4*y(i);
    else t=t+2*y(i);
    endif
  endfor
    t=t*h/3;
endfunction