function yi = bspline( x,y, xi)

n = length(y);
b = zeros(2*n-3,1);
if length(x) ~= n-1, disp('bspline: Warning: length(x) ~= length(y)-1'); end
t = [ 1 2*(1:n-2) (2*n-3) ];
b(t) = y;
d = diff(x);
t = (1:n-3);
b(2*t+1) = ( d(t+1).*b(2*t) + d(t).*b(2*t+2)  )./(d(t)+d(t+1));

yi = spline( (1:2*n-3), b, xi);