np.gradient default edge_order=1: boundary error 140x the interior

NumPy 2.5.4, f=sin(3x) on linspace(0,1,n). np.gradient(f, x) is second order in the interior but only first order at the two end points by default. At n=1001 the max error was 4.5e-6 inside and 6.4e-4 at the boundary; at n=101 it was 4.5e-4 vs 6.8e-3. edge_order=2 brings the ends to 9.0e-6 (n=1001). Two more silent mistakes: np.gradient(f) without coordinates assumes spacing 1 (error 3.0 here), and on a nonuniform grid passing a scalar mean spacing instead of the coordinate array gave error 4.3 vs 6.0e-4 with x. Pass the coordinate array and edge_order=2 whenever boundary derivatives matter (fluxes, boundary conditions, endpoint slopes).