function [y, yp, ypp] = bsplines(x)

% [y, yp, ypp] = bsplines(x)
%
% Computes the values and the first and second derivatives of the 13 cubic 
% B-splines with breakpoints 0, 1, 2, ..., 10.   
% x must be in the interval [0, 10].
%
% y:    a column vector with the values of the 13 cubic B-splines evaluated
%       at x.
% yp:   a column vector with the first derivatives of the 13 cubic B-splines
%       B-splines evaluated at x.
% ypp:  a column vector with the second derivatives of the 13 cubic 
%       B-splines evaluated at x.

if x < 0 || x > 10
   error('x must be in [0, 10]');
end;

% knots
a = [ 0, 0, 0, 0:10, 10, 10, 10 ];

% 0 <= x < 10 and a(i) <= x < a(i+1) or x = 10 and i is 14.
i = floor(x) + 4;

% y1, y2, y are vectors of linear, quadratic, and cubic B-splines at x.
if i == 14, 
    y1 = [ zeros(10,1); 1 ];
    y2 = [ zeros(11,1); 1 ];
    y = [ zeros(12,1); 1 ];
else
    w1 = (x - a(i)) / (a(i+1) - a(i));
    G1 = [ 1 - w1; w1 ];
    w2 = [ (x - a(i-1)) / (a(i+1) - a(i-1) ) ; 
           (x - a(i)) / (a(i+2) - a(i)) ];
    G2 = [ diag(1-w2); 0, 0] + [ 0, 0; diag(w2) ];
    w3 = [ (x - a(i-2)) / (a(i+1) - a(i-2)) ;
           (x - a(i-1)) / (a(i+2) - a(i-1)) ;
           (x - a(i)) / (a(i+3) - a(i)) ];
    G3 = [ diag(1-w3); 0, 0, 0 ] + [ 0, 0, 0; diag(w3) ];
    y = zeros(13, 1);
    y(i-3:i) = G3*G2*G1; 
    y2 = zeros(12, 1);
    y2(i-3:i-1) = G2*G1;
    y1 = zeros(11, 1);
    y1(i-3:i-2) = G1;
end;

% Width of support of linear, quadratic, and cubic B-splines.
d1 = a(5:15)' - a(3:13)';  
d2 = a(5:16)' - a(2:13)';  
d3 = a(5:17)' - a(1,13)';  

% Derivatives of quadratic B-splines.
y2p = zeros(12,1);
y2p(1) = -2 * y1(1) / d1(1);
y2p(2:11) = 2 * ( y1(1:10) ./ d1(1:10) - y1(2:11) ./ d1(2:11) );
y2p(12) = 2 * y1(11) / d1(11);

% First derivatives of cubic B-splines.
yp = zeros(13,1);
yp(1) = -3 * y2(1) / d2(1);
yp(2:12) = 3 * ( y2(1:11) ./ d2(1:11) - y2(2:12) ./ d2(2:12) );
yp(13) = 3 * y2(12) / d2(12);

% Second derivatives of cubic B-splines.
ypp = zeros(13,1);
ypp(1) = -3 * y2p(1) / d2(1);
ypp(2:12) = 3 * ( y2p(1:11) ./ d2(1:11) - y2p(2:12) ./ d2(2:12) );
ypp(13) = 3 * y2p(12) / d2(12);