function  [int,err]=simpson(f,a,b,n)
%
% function  [int,err]=simpson(f,a,b,n)
%  Computes Simpson's Rule with error estimate. 
%  Assumes that n is divisible by 4.
%
       n4=n/4; n2=n/2;
       int4=trap(f,a,b,n4);
       [int2,err2]=update_trap(f,a,b,n2,int4);
       S2=int2+err2;
       [int0,err0]=update_trap(f,a,b,n,int2);
       int=int0+err0;
       err=(int-S2)/15;