BUG: Fixes #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 `second_order` keyword (Default: False) is added to specify whether first or second order calculations at the boundaries should be used.

Tests and documentation updated to reflect this change.