function [x,w,varargout] = simprule(n,a,b)
%
% Generate nodes and weights for
% Simpson's rule on [a,b].
% Note that x is a column vector
% and w is a row vector.  Note also
% that n is the number of points, not
% the number of intervals, which is n-1.
%
nw = n;
if nw == 2*floor(n/2) % make sure n is odd
   nw = nw+1;
end
x = linspace(a,b,nw).';
w = repmat((b-a)/(nw-1)/3,1,nw);
w(2:end-1) = 2*w(2:end-1);
w(2:2:end-1) = 2*w(2:2:end-1);
if nargout==3
    varargout{1} = nw;
end