function c = fdstencil(k,j)

% Compute stencil coefficients for finite difference approximation 
% of k'th order derivative on a uniform grid.  Print the stencil and
% dominant terms in the truncation error.
%
% j should be a column vector
% of indices of grid points, values u(x0 + j*h)
% are used, where x0 is an arbitrary grid point and h the mesh spacing.
% This routine returns a vector c of length n=length(j) and the
% k'th derivative is approximated by
%   1/h^k * [c(1)*u(x0 + j(1)*h) + ... + c(n)*u(x0 + j(n)*h)].
% Typically j(1) <= 0 <= j(n) and the values in j are monotonically
% increasing, but neither of these conditions is required.
% The routine fdcoeffF is used to compute the coefficients.
%
% Example:   fdstencil(2,-1:1);
%   determines the 2nd order centered approximation of the 2nd derivative.
%
% From  http://www.amath.washington.edu/~rjl/fdmbook/  (2007)

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

c = fdcoeffF(k,0,j);      % coefficients for k'th derivative 

% print out stencil:

disp(' ')
disp(sprintf('The derivative u^(%i) of u at x0 is approximated by',k))
disp(' ')
disp(sprintf('          1/h^%i * [',k))
for i=1:n-1
  if j(i) < 0
    disp(sprintf('                     %22.15e * u(x0%i*h) + ',c(i),j(i)))
  elseif j(i) == 0
    disp(sprintf('                     %22.15e * u(x0) + ',c(i)))
  else
    disp(sprintf('                     %22.15e * u(x0+%i*h) + ',c(i),j(i)))
  end
end

disp(sprintf('                     %22.15e * u(x0+%i*h) ]   ',c(n),j(n)))


% Determine dominant terms in truncation error and print out:

err0 = c*(j(:).^n) / factorial(n);
err1 = c*(j(:).^(n+1)) / factorial(n+1);
if (abs(err0)) < 1e-14,  err0 = 0; end   % for centered approximations, expect
if (abs(err1)) < 1e-14,  err1 = 0; end   % one of these to be exactly 0.
disp(' ')
disp('For smooth u,')
disp(sprintf('       Error = %g * h^%i*u^(%i)(x0) + %g * h^%i*u^(%i)(x0) + ...',err0,n-k,n,err1,n-k+1,n+1))

disp(' ')


% Modifications to obtain fractions for the coefficients c(i).

% print out stencil:

disp(' ')
disp(sprintf('The derivative u^(%i) of u at x0 is approximated by',k))
disp(' ')
disp(sprintf('          1/h^%i * [',k))
for i=1:n-1
  if j(i) < 0
    fprintf('%s * u(x0%i*h) + ',strtrim(rats(c(i))),j(i))
  elseif j(i) == 0
    fprintf('%s * u(x0) + ',strtrim(rats(c(i))))
  else
    fprintf('%s * u(x0+%i*h) + ',strtrim(rats(c(i))),j(i))
  end
end

disp(sprintf('%s * u(x0+%i*h) ]   ',strtrim(rats(c(n))),j(n)))

% Dominant terms in truncation error with a fraction as a coeffcient:

disp(' ')
disp('For smooth u,')
disp(sprintf('       Error = %s * h^%i*u^(%i)(x0) + %s * h^%i*u^(%i)(x0) + ...',strtrim(rats(err0)),n-k,n,strtrim(rats(err1)),n-k+1,n+1))