function a = divided_difference( x, y )
  % ----------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 7, 1999
  %
  % Usage: a = divided_difference(x, y)
  %
  % Returns:
  %	array containing divided differences (polynomial coefficients)
  %	
  % Parameters:
  %	x:	array of x-values
  %	y:	array of y-values
  %
  % This function computes the divided differences that can be used
  % to construct an interpolating polynomial for the given set of data.
  %
  % This function accepts two vectors and returns a vector whose elements
  % are the divided difference values.
  %
  % The value of the polynomial p (degree will be at most n) evaluated at r
  % is given by:
  %
  %	s = a(n+1);
  %	for i = n : -1 : 1
  %	    s = s * ( r - x(i) ) + a(i);
  %	end
  %
  % using the vector a() returned by this function.
  % ----------------------------------------------------------------------

  [nx, mx] = size( x );
  [ny, my] = size( y );

  if ( max( nx, mx ) != max( ny, my ) )
    fprintf( 'both input variables must be vectors of the same length\n' );
    return;
  elseif ( ( min( nx, mx ) != 1 ) || ( min( ny, my ) != 1 ) )
    fprintf( 'both input variables must be vectors, not matrices\n' );
    return;
  end

  n = max( nx, mx ) - 1;
  a = y;

  for i = 1 : n
    for j = n+1 : -1 : i+1
      a(j) = ( a(j) - a(j-1) ) / ( x(j) - x(j-i) );
    end
  end

end
