Simpson Rule issue - Approximation worsens as n increases

I'm trying to create a function to estimate function area via the 1/3 Simpson Rule, but my results are odd. When I set the number of subintervals (n) to 2 (the lowest possible number for Simpson's Rule), it gets near the real answer. However as I increase n, my approximation gets worse (which is the opposite of what should happen). Can anyone help find what's wrong with my code? Thanks!

function sum=CompositeSimpson(f,a,b,n)
if (mod(n,2)~=0 ) 
    disp('User! Give me an even number of subintervals!'); 
    area=NaN; return; 
end
sum=0;
h=(b-a)/n;
for i=1:(n/2)  
    a1=a+2*(i-1)*h;
    b1=a+2*i*h;
    sum=sum+((b1-a1)/2)/3*(f(a1)+4*f((b1-a1)/2)+f(b1));

end
end

My testing is using f=@(x) x^3-x^2+x, a=0, b=1. Here are the n values I used and their results:

2 - 0.4167

4 - 0.2812

6 - 0.2377

8 - 0.2148

10 - 0.2007

16 - 0.1785

20 - 0.1709

40 - 0.1552

Answer 1

[Jim Riggs]

I have not tried to run your code, but at first glance I see a typo in the code.

sum=sum+((b1-a1)/2)/3*(f(a1)+4*f((b1-a1)/2)+f(b1));

In this line, the middle function evaluation is supposed to be at the mid point between a and b, so the expression should be:

sum=sum+((b1-a1)/2)/3*(f(a1)+4*f((b1+a1)/2)+f(b1));

This can also be simplified to

sum=sum+(b1-a1)/6*(f(a1) + 4*f((b1+a1)/2) + f(b1));

I have not yet checked your loop indexing to verify that it is correct.

[Jim Riggs]

After further review, it looks like your code should work with the above change.

[John D'Errico]

+1. Yes, it should work then. With only the additional caveat that using the variable name sum is just a friendly invitation to further bugs down the road, when you might actually need to use the function sum. Of course, then Andrew will ask the frenzied question about why the function sum does not seem to work for him. Naming variables to match the names of tremendously useful functions is a bad idea.

Others might note that the use of a variable named i MIGHT be another subtle problem in the future, in case you ever need to work with complex data. I tend to use i myself though, a holdover from the dark days of Fortran, so I can't complain too much. Even so, the use of i can be dangerous to those who are just learning MATLAB.

[Jim Riggs]

Thanks John, you raise some good points. (I also am an old FORTRAN guy myself and am quite comfortable using i as a loop index. I appreciate your reminder.)