>>> x = np.array([0., 2., 3., 6., 8.])
>>> y = x ** 2
>>> np.gradient(y, x, edge_order=2)
array([ 0.,  4.,  6., 12., 16.])

2. A scalar, representing the fixed sample distance:

>>> dx = 2
>>> x = np.array([0., 2., 4., 6., 8.])
>>> y = x ** 2
>>> np.gradient(y, dx, edge_order=2)
array([ 0.,  4.,  8., 12., 16.])

It's possible to provide different data for spacing along each dimension.
The number of arguments must match the number of dimensions in the input
data.

>>> dx = 2
>>> dy = 3
>>> x = np.arange(0, 6, dx)
>>> y = np.arange(0, 9, dy)
>>> xs, ys = np.meshgrid(x, y)
>>> zs = xs + 2 * ys
>>> np.gradient(zs, dy, dx)  # Passing two scalars
(array([[2., 2., 2.],
        [2., 2., 2.],
        [2., 2., 2.]]),
 array([[1., 1., 1.],
        [1., 1., 1.],
        [1., 1., 1.]]))

Mixing scalars and arrays is also allowed:

>>> np.gradient(zs, y, dx)  # Passing one array and one scalar
(array([[2., 2., 2.],
        [2., 2., 2.],
        [2., 2., 2.]]),
 array([[1., 1., 1.],
        [1., 1., 1.],
        [1., 1., 1.]]))

Notes
-----
Assuming that :math:`f\in C^{3}` (i.e., :math:`f` has at least 3 continuous
derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
minimize the "consistency error" :math:`\eta_{i}` between the true gradient
and its estimate from a linear combination of the neighboring grid-points:

.. math::

    \eta_{i} = f_{i}^{\left(1\right)} -
                \left[ \alpha f\left(x_{i}\right) +
                        \beta f\left(x_{i} + h_{d}\right) +
                        \gamma f\left(x_{i}-h_{s}\right)
                \right]

By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
with their Taylor series expansion, this translates into solving
the following the linear system:

.. math::

    \left\{
        \begin{array}{r}
            \alpha+\beta+\gamma=0 \\
            \beta h_{d}-\gamma h_{s}=1 \\
            \beta h_{d}^{2}+\gamma h_{s}^{2}=0
        \end{array}
    \right.

The resulting approximation of :math:`f_{i}^{(1)}` is the following:

.. math::

    \hat f_{i}^{(1)} =
        \frac{
            h_{s}^{2}f\left(x_{i} + h_{d}\right)
            + \left(h_{d}^{2} - h_{s}^{2}\right)f\left(x_{i}\right)
            - h_{d}^{2}f\left(x_{i}-h_{s}\right)}
            { h_{s}h_{d}\left(h_{d} + h_{s}\right)}
        + \mathcal{O}\left(\frac{h_{d}h_{s}^{2}
                            + h_{s}h_{d}^{2}}{h_{d}
                            + h_{s}}\right)

It is worth noting that if :math:`h_{s}=h_{d}`
(i.e., data are evenly spaced)
we find the standard second order approximation:

.. math::

    \hat f_{i}^{(1)}=
        \frac{f\left(x_{i+1}\right) - f\left(x_{i-1}\right)}{2h}
        + \mathcal{O}\left(h^{2}\right)

With a similar procedure the forward/backward approximations used for
boundaries can be derived.

References
----------
.. [1]  Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
        (Texts in Applied Mathematics). New York: Springer.
.. [2]  Durran D. R. (1999) Numerical Methods for Wave Equations
        in Geophysical Fluid Dynamics. New York: Springer.
.. [3]  Fornberg B. (1988) Generation of Finite Difference Formulas on
        Arbitrarily Spaced Grids,
        Mathematics of Computation 51, no. 184 : 699-706.
        `PDF <https://www.ams.org/journals/mcom/1988-51-184/
        S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
Nr