function jumps=fnjmp(f,x)
%
%        jumps=fnjmp(f,x)
%
%  returns the matrix  f(x+) - f(x-)  of jumps at the points in  x  in the
% function specified by  f .

% copyright 1994 Carl de Boor
%  cb : 20 fe 1994

%  If necessary, sort  x :
xs = x(:)'; lx = length(xs); tosort = 0;
if (~isempty(find(diff(xs)<= 0))),
   tosort = 1; [xs,indx] = sort(xs);
end

if (f(1)==11) 
   [t,a,n,k,d] = spbrk(f);
   jumps = zeros(d,lx);
   indexr = sorted(t,xs); indexl = n+k - sorted(-t,-xs); 
   indexl = indexl(lx:-1:1);
   index = find(indexr-indexl>=k);
   if ~isempty(index)
      a = [zeros(d,1) a zeros(d,1)];
      jumps(:,index) = a(:,indexr(index)+2-k) - a(:,indexl(index)+1);
   end
else
   [breaks,coefs,l,k,d] = ppbrk(f);
   jumps = zeros(d,lx);
   if l==1, return, end
   index = find(xs-breaks(sorted(breaks(2:l),xs)+1)==0);
   if ~isempty(index)
      jumps(index) = ppual(f,xs(index)) - ppual(f,xs(index),1);
   end
end

if tosort>0, [ignored,indx] = sort(indx); jumps = jumps(:,indx); end