  function alpha = BVP_Collocation(a,ua,b,ub,q,r,n)
% a<b and ua and ub are scalars.
% q and r are handles to functions q(x) and r(x) defined on [a,b].
% n is a positive integer >=2.
% alpha is a column (n+2)-vector with the property that if
%
%           utilde(x) = alpha(1)B_0(x) + ... + alpha(n+2)*B_{n+1}(x)
% 
% then 
%           utilde(a) = ua
%          -utilde''(z(i)) + q(z(i))u(z(i)) = r(z(i))     i=1:n
%           utilde(b) = ub
%
% where z = linspace(a,b,n) and
% B_{k}(z) is the B-spline Bstar((z-xk)/h) where xk = a+(k-1)h
% and h = (b-a)/(n-1).

% Set up an (n+2)-by-(n+2) linear system T*alpha = rhs...

z = linspace(a,b,n)';
h = (b-a)/(n-1);
rhs = zeros(n+2,1);
T = zeros(n+2,n+2);
% First equation deals with th eleft end condition...
%        utilde(a) = ua
% i.e.,
%        alpha(1)*B_{0}(a) + ... + alpha(n+2)*B_{n+1)(a) = ua
rhs(1) = ua;
for k=0:n+1
    xk = a+(k-1)*h;
    T(1,k+1) = derBstar((a-xk)/h);
end

% The next n equations are
%        -utilde''(z(i)) + q(z(i))utilde(z(i)) = r(z(i))     i=1:n
% i.e.,
%  sum_{j=0}^{n+1} ( -B"_{j}(z(i)) + q(z(i))B_{j}(z(i)) ) alpha(j+1) = r(z(i))
for i=2:n+1
    zi = z(i-1);
    rhs(i) = r(zi);
    for k=0:n+1
        xk = a+(k-1)*h;
        [Bk,dBk,ddBk] = derBstar((zi-xk)/h);
        T(i,k+1) =  -ddBk/h^2 + q(zi)*Bk;
    end
end

% Last equation deals with the right end condition...
%        utilde(b) = ub
% i.e.,
%        alpha(1)*B_{0}(b) + ... + alpha(n+2)*B_{n+1)(b) = ub
rhs(n+2) = ub;
for k=0:n+1
    xk = a+(k-1)*h;
    T(n+2,k+1) = derBstar((b-xk)/h);
end
format short
alpha = T\rhs;