function b = bspline(xi,x) 
% BSPLINE: evaluation of a B-spline  
% b = bspline(xi,x) 
% values b(j) of a B-spline with knots xi at x(j)

b = 0*x; 
if length(xi) == 2
   % constant B-spline
   b( x>=xi(1) & x<xi(2) ) = 1;
else
   % recursion, discarding terms with zero denominator
   d = xi(end-1)-xi(1);
   if d > 0
      b = b + (x-xi(1)).*bspline(xi(1:end-1),x)/d;
   end
   d = xi(end)-xi(2);
   if d > 0
      b = b + (xi(end)-x).*bspline(xi(2:end),x)/d;
   end
end

