# Takes a function, start point, end point, and number of intervals.
# Returns a numerical approximation of the integral using Simpsons Rule.
sub integrate_simpson {
  my ($func, $start, $end, $n) = @_;

  my $width = ($end - $start)/$n;
  my $mid_sum = sum(map {$func->($start + ($_ - 0.5)*$width)} 1..$n);
  my $int_sum = sum(map {$func->($start + $_ *$width)} 1..$n-1);
  my $first = $func->($start);
  my $last = $func->($end);

  ($first + $last + 4*$mid_sum + 2*$int_sum) * $width / 6;
}

sub sum {
  my $sum = shift;
  $sum += shift while @_;
  return $sum;
}

# And to test
print integrate_simpson(sub {$_[0]**3}, 0, 1, 9)