function q = line_integral ( n )

%*****************************************************************************80
%
%% LINE_INTEGRAL...
%
  p1 = 0;
  p2 = 3;
  svec = rand ( n, 1 );
  x = p1(1) + svec * ( p2(1) - p1(1) );
  fvec = x .* x .* (4-x) + 1;
  q = norm ( p2 - p1 ) * sum ( fvec ) / n;

  return
end
