function T=simpson(f_name, a, b, n)
% f_name为要求的定函数y=f(x)所在的程序文件名
% a为积分下限
% b为积分上限
% n为积分区间[a,b]划分成小区间的等份数

h = (b-a) / n;
x = a + (0:n) * h;
f = feval(f_name, x);
N = length(f) - 1;
if N == 1
    fprintf('Data has only one interval')
    return;
end
if N == 2
    T = h / 3 * (f(1) + 4 * f(2) + f(3));
    return;
end
if N == 3
    T = 3 / 8 * h * (f(1) + 3 * f(2) + 3 * f(3) + f(4));
    return;
end
T = 0;
if 2 * floor(N / 2) == N
    T = h / 3 * (2 * f(N-2) + 2 * f(N-1) + 4 * f(N) + f(N+1));
    m = N - 3;
else
    m = N;
end
    T = T + (h/3) * (f(1) + 4 * sum(f(2:2:m)) + 2 * f(m+1));
if m > 2
    T = T + (h/3) * 2 * sum(f(3:2:m));
end

function y=simpson_f(x)
y = sin(x)
end

T = simpson('simpson_f', 0, pi, 4)

%=> y = 0    0.7071    1.0000    0.7071    0.0000
%=> T = 2.0046