I am trying to implement Simpson's rule, but for one reason or another, this works only for constant functions, but for every other kind of function, my code doesn't give a decent approximation. The code itself is rather simple, but I just don't see where I am mistaken. First, I built Simpson's rule as a function:
function [summe] = SimpsonQuad(x, f, n)
summe = 0;
for j = 1:n
summe = summe + ((x(j+2) - x(j))/6) * (f(x(j)) + 4*f(x(j+1)) + f(x(j+2)));
end
and then, I wrote a script to call it:
int = [0,3];
f = inline('x');
n = 10;
x = [];
x(1) = int(1);
for j = 1:2*n
x(j+1) = int(1) + (j*(int(2)-int(1))/(2*n));
end
SimpsonQuad(x, f, n);
As you can see, I calculate the sampling points $x_j$ in the script, and this works the way it should. But the result isn't decent. For example, for $n = 2$, I receive
0 0.750000000000000 1.500000000000000 2.250000000000000 3.000000000000000
ans =
3.375000000000000
For greather $n$, the result becomes even less, which doesn't make any sense at all.
I'd appreciate any help!
