function [sum] = simp(alpha,n)

%  Computes an approximation to the integral from 0 to 1 of
%  cos(alpha x^2) dx  using Simpson's rule with n points.

dx = 1/n;

sum = 0;
for i=0:n,              % Loop over nodes.
  x = i*dx;
  cosaxsq = cos(alpha*x^2);  % Evaluate integrand at x.
  if i==0 | i==n, 
    sum = sum + cosaxsq/6;   % Add 1/6 * Value of integrand at endpoints.
  else
    sum = sum + cosaxsq/3;   % Add 1/3 * Value of integrand at interior points.
  end;

  if i < n,
    xmid = x + .5*dx;
    cosaxmsq = cos(alpha*xmid^2);
    sum = sum + 2*cosaxmsq/3;     % Add 2/3 * Value of integrand at midpoints.
  end;
end;

sum = sum*dx;        %  Multiply final result by dx.