% plain and simple Simpson's rule

function Sf = simp(a,b,f,fa,fb)
global pts

c = (a+b)/2;
fc = feval(f,c);
pts = [pts; c];

Sf = (b-a)/6*(fa+4*fc+fb);