function [ df, A, B ] = vcompact4( tau, f, hin )
%VCOMPACT4 fourth-order compact finite difference derivative
%   Given a mesh of distinct points tau and a vector of function values on
%   tau use a fourth-order variable mesh compact finite difference formula 
%   compute df(k) = f'(tau(k)).  Except for the endpoints, this is the
%   same procedure as is used in SPLINE to compute derivatives making
%   the spline twice continuously differentiable.
%
%   [df, A, B ] = vcompact4( tau, f, h )
%   df = vector of derivative values
%   tau = vector of distinct nodes
%   f  = vector of function values
%   h (optional) vector of mesh widths or string 'cos' for Chebyshev-
%     Lobatto points
%
%   Optional output the matrices A and B for inspection.
%
%   RMC January 2011
%   
n = length( tau );
if nargin > 2,
    if hin == 'cos',
        % Accurate computation of widths, in the
        % common case of Chebyshev nodes cos(pi*(0:n-1)/(n-1)).
        m = floor( (n-1)/2 );
        s = sin( pi*(0:m)/(n-1) );
        % relative accuracy of computation of sine is possible for small
        % indexes but not near pi; use the Weidemann-Reddy flip trick
        % 
        if rem(n,2)==1,
            s = [s, fliplr(s(1:end-1))];
        else
            s = [s, fliplr(s)];
        end
        h = -s(2)*( s(1:n-1)+s(2:n) )/(1+cos(pi/(n-1)));
    else
        h = hin;
    end
else
    h = diff( tau ); 
end
A = sparse( n, n );
B = sparse( n, n );

% Special formula at the left end
A(1,1) = 1.0/(h(2)+h(1))/(h(3)+h(2)+h(1));
A(1,2) = 1.0/h(2)/(h(3)+h(2));
% Unless carefully written, these are quite susceptible to rounding errors.
B(1,1) = (4*h(1)^2+6*h(1)*h(2)+3*h(1)*h(3)+2*h(2)^2+2*h(2)*h(3)) ...
         /(h(1)+h(2))^2/(h(1)+h(2)+h(3))^2/h(1);
B(1,2) = 1/h(1)*((-2*h(2)+h(1))*h(3)+2*h(2)*(-h(2)+h(1))) ...
                 /h(2)^2/(h(2)+h(3))^2;
B(1,3) = -h(1)^2/(h(2)+h(1))^2/(h(2))^2/h(3);
B(1,4) =  h(1)^2/(h(3)+h(2)+h(1))^2/(h(3)+h(2))^2/h(3);

% Symmetric compact finite difference formula on interior nodes
for k=2:n-1,
    A(k,k-1) = 1/(h(k-1)+h(k))^2;
    A(k,k)   = 1/(h(k))^2;
    A(k,k+1) =  h(k-1)^2/(h(k-1)+h(k))^2/(h(k))^2;
    B(k,k-1) = 2*(2*h(k-1)+h(k))/h(k-1)/(h(k-1)+h(k))^3;
    B(k,k)   = -2*(-h(k-1)+h(k))/h(k-1)/h(k)^3;
    B(k,k+1) = -2*h(k-1)^2*(h(k-1)+2*h(k))/(h(k-1)+h(k))^3/h(k)^3;
end


% Special formula at the right end
A(n,n-1) = 1/(h(n-2))/(h(n-2)+h(n-3));
A(n,n)   = 1/(h(n-1)+h(n-2))/(h(n-3)+h(n-2)+h(n-1) );
% As above, unless care is used, these are quite susceptible to rounding errors.
B(n,n-3) = -h(n-1)^2/(h(n-3)+h(n-2)+h(n-1))^2/(h(n-3)+h(n-2))^2/h(n-3);
B(n,n-2) = h(n-1)^2/(h(n-2)+h(n-1))^2/h(n-2)^2/h(n-3);
B(n,n-1) = 1/h(n-1)*((2*h(n-2)-h(n-1))*h(n-3)+2*h(n-2)*(h(n-2)-h(n-1))) ...
           /h(n-2)^2/(h(n-3)+h(n-2))^2;
B(n,n) = -(4*h(n-1)^2+(6*h(n-2)+3*h(n-3))*h(n-1)+2*h(n-2)*(h(n-3)+h(n-2))) ...
           /h(n-1)/(h(n-2)+h(n-1))^2/(h(n-3)+h(n-2)+h(n-1))^2;

% Solve the tridiagonal system for the unknown derivatives
df = -A\(B*f(:));
