Simpson's Rule

So I have to write a script to evaluate an integral using Simpson's Rule, including odd values of N. Right now I am not getting the correct answer even for even values of N. Here is what I have:

% Ask for user input
% Lower bound (a)
a = input('What is your lower bound (a)?')

% Upper bound (b)
b = input('What is your upper bound (b)?')

% Subintervals
N = input('How many subintervals (N)?')

% Defining function
f = @(x) (2*x)/((x.^2)+1)

% Finding h
h=(b-a)/N; 

% Finding the values of x for each interval
x=linspace(a,b,N);

% Calculating the integral
for i = 1:N-1
I(i)= (h/3)*(f(x(i))+(4*f((x(i)+x(i+1))/2))+f(x(i+1)));
end

answer1 = sum(I)

I'm really not sure where I'm going wrong. I have written it out on paper many times and still cannot find my error.

[Muhammad Mahboob Khalil]

Here is the fit script for Simpson's Rule:

% MATLAB code for syms function that creates a variable

% dynamically and automatically assigns

% to a MATLAB variable with the same name

syms x

% Lower Limit

a = 4;

% Upper Limit

b = 5.2;

% Number of Segments

n = 6;

% Declare the function

f1 = log(x);

% inline creates a function of string containing in f1

f = inline(f1);

% h is the segment size

h = (b - a)/n;

% X stores the summation of first and last segment

X = f(a)+f(b);

% variables Odd and Even to store

% summation of odd and even

% terms respectively

Odd = 0;

Even = 0;

for i = 1:2:n-1

xi=a+(i*h);

Odd=Odd+f(xi);

end

for i = 2:2:n-2

xi=a+(i*h);

Even=Even+f(xi);

end

% Formula to calculate numerical integration

% using Simpsons 1/3 Rule

I = (h/3)*(X+4*Odd+2*Even);

disp('The approximation of above integral is: ');

disp(I);

Answer 1

[Walter Roberson]

It looks to me as if your result is twice as large as it should be.

Notice in https://www.chegg.com/homework-help/definitions/simpsons-rule-29 that they give two forms. You are using the 1 4 1 form, for which the multiplier should be h/6 . The multiplier of h/3 is for the 1 4 2 4 2 ... 4 1 form

Answer 2

[Edoardo Bertolotti]

I am not sure about the version with

I(i)= (h/3)* [...]

but if you use

I(i)= (h/6) [...]

x=linspace(a,b,N+1);

for i = 1:N

it should work fine.