Not entirely clear to me between what limits you are integrating, nor what your loops are doing! Perhaps the following might help:
% Integrate exp(-x^2) from x = 0 to x = 3 using Simpson's 3/8 rule
% I = 3*h/8*(f(a))+sum(f(i)+f(i+1),i=2:3:n-2)+sum(f(i),i=3:3:n-1)+f(b))
f = @(x) exp(-x.^2); % function
a = 0; % lower limit
b = 3; % upper limit
n = 333; % number of panels (multiple of 3)
h = (b - a)/n; % panel size
x = a:h:b; % vector of x values
s3 = 0;
for j = 2:3:n-2
s3 = s3 + f(x(j)) + f(x(j+1));
end
s2 = 0;
for j = 3:3:n-1
s2 = s2 + f(x(j));
end
I = 3*h/8*(f(a) + 3*s3 + 2*s2 + f(b));
disp('Approximate integral of exp(-x^2) from x=0 to x=3 using')
disp(['Simpson''s 3/8 rule with ' int2str(n) ' panels is: ' num2str(I,4)])
Modify to suit exactly what you want.
