% Integrate the function in f.m using simpson's rule


N = 9;  %Number of points in integration
x1 = 0; % Lower limit of integration
xN = 5; % Upper limit of integration
delta = (xN-x1)/(N-1); 
integral = 0;
for n=1:1:N
	xn = x1 + (n-1)*delta;

	if n==1 | n==N
		integral = integral + f(xn);
	elseif mod(n,2)==0   %If the remainder of n/2 is zero (n is even)
		integral = integral + 4*f(xn);
	else
		integral = integral + 2*f(xn);
	end
   integral = (delta/3)*integral;
end
fprintf('The integral of f.m between %g and %g is %g \n',x1,xN,integral);