simpson.m

%
% Use Matlab's symbolic features to derive Simpson's 1/3 rule
%

% define matrix and vector...
%
%   assume data pairs of: (0, y1), (h, y2), (2h, y3)
%
%   fit to: f(x) = a_0 + a_1*x + a_2*x^2
%

A = sym('[1, 0, 0; 1, h, h^2; 1, 2*h, 4*h^2]')
b = sym('[y1; y2; y3]')

% solve system of equations

a = linsolve(A,b)

% form a_0 + a_1*x + a_2*x^2

f = symadd(symadd(sym(a,1,1),symmul(sym(a,2,1),'x')),symmul(sym(a,3,1),'x^2'))

% perform integration

i = int(f,0,'2*h')

% results
%
% A =  [ 1,   0,     0]
%      [ 1,   h,   h^2]
%      [ 1, 2*h, 4*h^2]
%
% b =  [ y1]
%      [ y2]
%      [ y3]
%
% a =  [                   y1]
%      [-1/2*(3*y1+y3-4*y2)/h]
%      [ 1/2*(y1+y3-2*y2)/h^2]
%
% f =  y1-1/2*(3*y1+y3-4*y2)/h*x+1/2*(y1+y3-2*y2)/h^2*x^2
%
% i =  1/3*h*(y1+4*y2+y3)   as we expect!
%