function I=simp4_1(Fun,a,b,n);
%%   function I=simp4_1(Fun,a,b,n);
%%   
%%   COMPOSITE SIMPSON'S RULE FOR INTEGRATION 
%%
%%   INPUT:
%%
%%    Fun - the function data which may be entered EITHER as
%%           a string, enclosed in ' ', with the name the m-file of the 
%%           function ,  
%%           OR as a row vector of data [f(x0),f(x1),...,f(xn)]
%%     a  - lower endpoint of integration
%%     b  - upper endpoint of integration
%%     n  - an even integer BIGGER THAN 2 (x_n is the endpoint b, altogether
%%                                         n+1 points counting  endpoints)
%%   OUTPUT:
%%     I - the approximate value of the integral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h=(b-a)/n;

if isstr(Fun)==1;  % test whether from a function file or from given data
   points=a+h*[0:n];
   Fvalues=feval(Fun,points);
end;

I=Fvalues(1)+Fvalues(n+1)+4*sum(Fvalues(2:2:n))+2*sum(Fvalues(3:2:n-1));
I=h*I/3;
