function c = fdcoeffV(k,xbar,x)

% Compute coefficients for finite difference approximation for the
% derivative of order k at xbar based on grid values at points in x.
%
% WARNING: This approach is numerically unstable for large values of n since
% the Vandermonde matrix is poorly conditioned.  Use fdcoeffF.m instead,
% which is based on Fornberg's method.
%
% This function returns a row vector c of dimension 1 by n, where n=length(x),
% containing coefficients to approximate u^{(k)}(xbar), 
% the k'th derivative of u evaluated at xbar,  based on n values
% of u at x(1), x(2), ... x(n).  
%
% If U is a column vector containing u(x) at these n points, then 
% c*U will give the approximation to u^{(k)}(xbar).
%
% Note for k=0 this can be used to evaluate the interpolating polynomial 
% itself.
%
% Requires length(x) > k.  
% Usually the elements x(i) are monotonically increasing
% and x(1) <= xbar <= x(n), but neither condition is required.
% The x values need not be equally spaced but must be distinct.  
%
% From  http://www.amath.washington.edu/~rjl/fdmbook/  (2007)

n = length(x);
if k >= n 
   error('***  length(x) must be larger than k')
   end

A = ones(n,n);
xrow = (x(:)-xbar)';  % displacements x-xbar as a row vector.
for i=2:n
   A(i,:) = (xrow .^ (i-1)) ./ factorial(i-1);
   end

b = zeros(n,1);    % b is right hand side,
b(k+1) = 1;        % so k'th derivative term remains

c = A\b;           % solve n by n system for coefficients
c = c';            % row vector