simp <- function ( f, a, b, m = 100 )

#*****************************************************************************80
#
## simp estimates an integral using Simpson's rule.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    04 March 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017
#    ISBN13: 978-1-4987-2363-3.
#
{
  x.ends <- seq ( a, b, length.out = m + 1 )
  y.ends <- f ( x.ends )

  x.mids <- ( x.ends[2:(m+1)] - x.ends[1:m] ) / 2.0 + x.ends[1:m]
  y.mids <- f ( x.mids )

  area = sum ( y.ends[2:(m+1)] + 4.0 * y.mids[1:m] + y.ends[1:m] )
  area <- area * abs ( b - a ) / 6.0 / m

  return ( area )
}