function spline = bkfspmak(x,j,Theta)

%
%      spline = bkfspmak(x,j,Theta)
%
%  make the (Birkhoff) spline which represents the j-th derivative at x of the
%  error in Hf the Hermite interpolant to f at Theta a multiset of n points, as:
%
%       D^j f(x) - D^j Hf(x) = integral spline(t)*D^n f(t) dt
%  
%  example:  spline = bkfspmak(0.1,2,[0 0 0 1 1 1])
%
%  See also:  bkfplt (which plots such a spline)
%
%  S. Waldron October 1994
%

n=length(Theta);   % number of points in Theta
r=n-length(find(Theta-x*ones(size(Theta))));  % r the multiplicity of x in Theta

if j>=n, error('j must be less than the number of points in Theta'); end

if j<r, spline = spmak(sort([Theta,x]),0); return; end
% the error is zero, and so spline is the zero spline

% the spline is to be expressed as: sum_{i=1}^{j+1-r} B(i)a(i)  
% where t is the knot sequence
t=sort([Theta,x*ones(1,j+1-r)]);
% and B(i) are the n-th order B-splines for t
% the spline satisfies j+1-r conditions, namely that the jump in the n-j-1
% derivative at x is (-1)^{n-1-j} and the spline does not contain the 
% truncated powers (.-x)_+^{n-j},...,(.-x)_+^{n-r-1}, i.e., has continuous
% n-j,...,n-r-1 derivatives at x. We write these (linearly independent)
% conditions as a linear system: C a = b, with a the vector of coefficients for
% the spline, solve for a, and make up the spline

% make the B-splines
bsplines = spmak(t,eye(j+1-r,j+1-r));
% compute C
bsplines=fnder(bsplines,n-j-1);
for i=1:j+1-r, C(i,:)=fnjmp(bsplines,x)'; bsplines=fnder(bsplines); end 
% make the corresponding b
b=zeros(j+1-r,1); b(1,1)=(-1)^(n-j-1);
% find a
a=C\b;
% make the spline
spline=spmak(t,a');

% hope for the best