% OTHERSIMP  Compute coefficients in second rule in 7.2 exercise #21.

% build the Lagrange polynomials as functions:
x = [0.25 0.5 0.75]; 
l0 = @(y) (y-x(2)) .* (y-x(3)) / ((x(1)-x(2)) * (x(1)-x(3)));
l1 = @(y) (y-x(1)) .* (y-x(3)) / ((x(2)-x(1)) * (x(2)-x(3)));
l2 = @(y) (y-x(1)) .* (y-x(2)) / ((x(3)-x(1)) * (x(3)-x(2)));

% plot them to check they are the right functions:
xf=0:0.001:1;
plot(x,ones(size(x)),'o',x,zeros(size(x)),'o'), hold on
plot(xf,l0(xf),xf,l1(xf),xf,l2(xf)), hold off, grid on
legend('l_0(x)','l_1(x)','l_2(x)')

% now apply Simpson's rule to *exactly* do the integrals which
%   determine the Newton-Cotes coefficients:
format rat
alpha = (1/6) * (l0(0) + 4*l0(0.5) + l0(1))
beta  = (1/6) * (l1(0) + 4*l1(0.5) + l1(1))
gamma = (1/6) * (l2(0) + 4*l2(0.5) + l2(1))
