function  int=trap(f,a,b,n)
%
%  function  int=trap(f,a,b,n)
%   Basic composite trapezoid rule for n points.
%
       h=(b-a)/n;
       fa=feval(f,a); fb=feval(f,b);
       int=0.5*h*(fa+fb);
       sum=0;
       for j=1:n-1
        xx=a+j*h;
        sum=sum+ feval(f,xx);
       end;
       int=int+h*sum;