function delta = differ_arbitrary ( M, N, x0, alpha )

%*****************************************************************************80
%
%% differ_arbitrary() approximates a derivative using arbitrary spaced data.
%  
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 September 2024
%
%  Author:
%
%    Original pseudocode by Bengt Fornberg.
%    This version by John Burkardt.
%
%  Reference:
%
%    Bengt Fornberg,
%    Generation of finite difference formulas on arbitrarily spaced grids,
%    Mathematics of Computation,
%    Volume 51, Number 184, October 1988, pages 699-706.
%
%  Input:
%
%    integer M: the order of the highest derivative to approximate.
%
%    integer N: there are N+1 grid point values.
%
%    real x0: the point at which the approximation is to hold.
%
%    real alpha(n+1): the x coordinates of the grid point values.
%    These should be distinct, but need not be evenly spaced.
%
%  Output:
%
%    real delta():
%
  delta = zeros ( min(N+1,M+1), N+1, N+1 );
%
%  O is an offset, typically either 1 or 0.
%
  O = 1;

  delta(0+O,0+O,0+O) = 1.0;

  c1 = 1.0;

  for n = 1 : N

    c2 = 1.0;

    for nu = 0 : n - 1

      c3 = alpha(n+O) - alpha(nu+O);
      c2 = c2 * c3;

      for m = 0 : min ( n, M )

        if ( m == 0 )
          delta(m+O,n+O,nu+O) = ...
            ( alpha(n+O) - x0 ) * delta(m+O,n-1+O,nu+O) / c3;
        else
          delta(m+O,n+O,nu+O) = ( ...
            ( alpha(n+O) - x0 ) * delta(m+O,n-1+O,nu+O) ...
            - m * delta(m-1+O,n-1+O,nu+O) ...
            ) / c3;
        end
      end
    end

    for m = 0 : min ( n, M )
      if ( m == 0 )
        delta(m+O,n+O,n+O) = c1 * ( ...
          - ( alpha(n-1+O) - x0 ) * delta(m+O,n-1+O,n-1+O) ...
          ) / c2;
      else
        delta(m+O,n+O,n+O) = c1 * ( ...
          m * delta(m-1+O,n-1+O,n-1+O) ...
          - ( alpha(n-1+O) - x0 ) * delta(m+O,n-1+O,n-1+O) ...
          ) / c2;
      end
    end

    c1 = c2;

  end

  return
end
  
