function s=simpson(f,a,b,n)
    if mod(n,2)~=0
       fprintf('simpson: n=%g is not even.\n',n);
       return;
    end;
    h=(b-a)/n;    
    so=0;
    for i=1:2:n-1
       so=so+f(a+h*i);
    end
    se=0;
    for i=2:2:n-2
       se=se+f(a+h*i);
    end
    s=h*(f(a)+f(b)+4*so+2*se)/3;
