Simpson's Rule

I coded Simpson's Rule, but I'm not sure if it's right.

f = @(x) exp(-x.^2);
true = integral(f, 0, 1);
%simpson's rule
n = 128;
k= n/2;
a = 0; b = 1;
h = (b-a)/n;
x = a + h;
sum1 = (h/3)*(f(b)+ f(a));
sum1 =  sum1 + (4*h/3)*f(x);
for j = 1:n-1
    x1 = a+(2*j)*h;
    x2 = a + (2*j+1)*h;
    sum1 = sum1 + (2*h/3)*f(x1)+(4*h/3)*f(x2);
end

error1 = abs(true-sum1);

fprintf(' Simpsons %d : %0.8f \n', n, sum1);
fprintf('Simpsons Error: %0.8f \n', error1);

[Geoff Hayes]

Chris - what makes you think that the algorithm has been coded incorrectly? Presumably you must have a set of test data that you will use to validate the above. What do you notice when you do so?

Answer 1

[Roger Stafford]

I think you made an error on the line

   for j = 1:n-1

It should be

   for j = 1:k-1

where k = n/2. As it stands now, the x2 reaches a value of a+(2*n-1)*h which is far beyond the range from 0 to 1.