function z = gaussian2(f,a,b)
% GAUSSIAN2  Use n=2 Gauss-Legendre rule to approximate integral.  Note
% that c_1 = c_2 = 1.
%
% Example:
%   >> gaussian2(@(x) sin(x),0,pi/4)
%   >> exact = 1 - sqrt(2)/2          % = 0.2928932188

sh = @(t) 0.5 * ((b-a) * t + a+b);  % does shift and scale: x = sh(t)

t1 = - sqrt(3)/3;   t2 = - t1;

z = 0.5 * (b-a) * ( feval(f,sh(t1)) + feval(f,sh(t2)) );
