BUG: Fixed #5184 gradient calculation behavior at boundaries

* Previous expected behavior was that the gradient is computed using central
  differences in the interior and first differences at the boundaries.

* gradient was updated in v1.9.0 so that second-order accurate calculations are
  done at the boundaries, but this breaks expected behavior with old code, so
  `edge_order` keyword (Default: 1) is added to specify whether first or second
  order calculations at the boundaries should be used.

* Since the second argument is *varargs, in order to maintain compatibility
  with old code and compatibility with python 2.6 & 2.7, **kwargs is used, and
  all kwargs that are not `edge_order` raise an error, listing the offending
  kwargs.

* Tests and documentation updated to reflect this change.