Sipmson's rule in Matlab

I'm trying to program the following Simpson's rule problem:

This is my attempt so far:

    function simp = mysimp(fName,x2,x1,m)
    h=(x2-x1)/(m);
    x = x1:h:x2;
    f = fName(x);
    simp = f(1) + f(end);
    for I = 2:2:m
        simp = simp + 4*(f(I));
    end
    for I = 3:2:m
        simp = simp + 2*(f(I));
    end
    simp = simp*(h/3);             
 end

This isn't returning the right answer for me. Can someone help me out?

Answer 1

[Geoff Hayes]

Ollie - I think that you almost have the correct implementation for Simpson's algorithm, but you may want to reconsider how you are determining h. Note that the question is telling you what h should be. So is your calculation giving you the same value? You may want to pass h into your function and calculate h from that. Like the following

 func = @(x) 2*x + 2*sin(x) - 4;
 mysimp(func, 3.1,3,0.1)

and your function mysimp is defined as

 function simp = mysimp(fName,x2,x1,h)

 m = floor((x2-x1)/h);

 % etc.

Note how we calculate m from the inputs and so have m subintervals. Now notice how you are creating the array of x values

 x = x1:h:x2;

A problem with this is that you are not guaranteed that the last value in this array will be x2 (it all depends on x1 and the step size h). Instead, you should use the linspace function as

 x = linspace(x1,x2,m);

to create an array of m+1 elements starting from x1 and ending at x2 (the step size between each value will be near enough to h). We then calculate the y values as

 y = fName(x);