function c = fdcoeffF(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.
%
% 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.
%
% This program should give the same results as fdcoeffV.m, but for large
% values of n is much more stable numerically.
%
% Based on the program "weights" in
% B. Fornberg, "Calculation of weights in finite difference formulas",
% SIAM Review 40 (1998), pp. 685-691.
%
% Note: Forberg's algorithm can be used to simultaneously compute the
% coefficients for derivatives of order 0, 1, ..., m where m <= n-1.
% This gives a coefficient matrix C(1:n,1:m) whose k'th column gives
% the coefficients for the k'th derivative.
%
% In this version we set m=k and only compute the coefficients for
% derivatives of order up to order k, and then return only the k'th column
% of the resulting C matrix (converted to a row vector).
% This routine is then compatible with fdcoeffV.
% It can be easily modified to return the whole array if desired.
%
% From http://www.amath.washington.edu/~rjl/fdmbook/ (2007)
n = length(x);
if k >= n
error('*** length(x) must be larger than k')
end
m = k; % change to m=n-1 if you want to compute coefficients for all
% possible derivatives. Then modify to output all of C.
c1 = 1;
c4 = x(1) - xbar;
C = zeros(n-1,m+1);
C(1,1) = 1;
for i=1:n-1
i1 = i+1;
mn = min(i,m);
c2 = 1;
c5 = c4;
c4 = x(i1) - xbar;
for j=0:i-1
j1 = j+1;
c3 = x(i1) - x(j1);
c2 = c2*c3;
if j==i-1
for s=mn:-1:1
s1 = s+1;
C(i1,s1) = c1*(s*C(i1-1,s1-1) - c5*C(i1-1,s1))/c2;
end
C(i1,1) = -c1*c5*C(i1-1,1)/c2;
end
for s=mn:-1:1
s1 = s+1;
C(j1,s1) = (c4*C(j1,s1) - s*C(j1,s1-1))/c3;
end
C(j1,1) = c4*C(j1,1)/c3;
end
c1 = c2;
end
c = C(:,end)'; % last column of c gives desired row vector