% example for hw problem with simplified assumptions
% p(t) q(t) w(t) are assumed to be constants
% basic functions assumed to be B-spline
% May contain errors or typos
% 
n = 5; % number of collocation point including boundary
a = 0.0; % left boundary
b = 1; % right boundary
alpha = 0; % left boundary
beta = 1; % right boundary
h = (b-a)/(n-3); % book equation (14)
t = [a:h:b]; % book equation (15)
p = 0; % should be function of t
q = 0; % should be function of t
w = 2; % should be function of t
%%
v = zeros(n,1);
A = zeros(n,n);
% set up linear system, see equation in hw problem, or book equaiton (12) (13)
%             A                      c      =    v
%|(Lv_1)(t_1) ... (Lv_n)(t_1)     | |c_1  |   | w(t_1)  |
%|(Lv_1)(t_2) ... (Lv_n)(t_2)     | |c_2  | = | w(t_2)  |
%| ...                 ...        | |...  |   | ...     |
%|(Lv_1)(t_n-2) ... (Lv_n)(t_n-2) | |c_n-2|   | w(t_n-2)|
%|( v_1)(a)   ... ( v_n)(a)       | |c_n-1|   | alpha   | (left boundary condition)
%|( v_1)(b)   ... ( v_n)(b)       | |c_n  |   | beta    | (right boundary condition)  

% i=1 to n-2, collocation points
for i=1:n-2
    v(i) = w;
    for j=1:n
        % v_j(t) = B((t-a)/h-j+2) book equation (17)
        % Lu = u''(t) + p u' + q u, where u is v_j, t is t_i
        A(i,j) = B((t(i)-a)/h-j+2,2) + p*B((t(i)-a)/h-j+2,1) + q*B((t(i)-a)/h-j+2,0);
    end 
end
% i=n-1, left boundary
for j=1:n
    A(n-1,j) = B((a-a)/h-j+2,0);
end
v(n-1) = alpha;

% i=n, right boundary
for j=1:n
    A(n,j) = B((b-a)/h-j+2,0);
end
v(n) = beta;

% solve Ac = v
c = A\v;

%% residual
tt = linspace(a,b,2*n-5);
Lu = zeros(size(tt));
for i = 1:length(tt)
    %for each tt, calculate Lu = c_1(Lv_1)(tt) +...+ c_n(Lv_n)(tt)
    for j = 1:n
        Lu(i) = Lu(i) + c(j)*(B((tt(i)-a)/h-j+2,2) + p*B((tt(i)-a)/h-j+2,1) + q*B((tt(i)-a)/h-j+2,0));
    end
end
% evaluate residual = Lu-w 
res = Lu-w;
plot(tt,res,'-o')
% plot(xx,2*sin(xx),'-o',xx,yy)

%% B-spline
% d = 0,1,2 for derivitives 
function [y] = B(x,d)
    % book equation (16)
    % @ for anonymous function https://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html
    f1 = @(t)  (t+2)^3/6;
    f2 = @(t)  (1+3*(t+1)+3*(t+1)^2-3*(t+1)^3)./6;
    % first derivatives
    d1f1 = @(t)  (t+2)^2/2;
    d1f2 = @(t)  (3+6*(t+1)-9*(t+1)^2)./6;
    % second derivatives
    d2f1 = @(t)  (t+2);
    d2f2 = @(t)  (6-18*(t+1))./6;

    switch d
        case 0
            for i=1:length(x)
                if abs(x(i))>2
                    y(i) = 0;
                elseif x(i)<=-1 && x(i)>=-2
                    y(i) = f1(x(i));
                elseif x(i)<=0 && x(i)>=-1
                    y(i) = f2(x(i));
                elseif x(i)>=0 && x(i)<=1
                    y(i) = f2(-x(i));
                elseif x(i)<=2 && x(i)>=1
                    y(i) = f1(-x(i));
                end
            end
        case 1
            for i=1:length(x)
                if abs(x(i))>2
                    y(i) = 0;
                elseif x(i)<=-1 && x(i)>=-2
                    y(i) = d1f1(x(i));
                elseif x(i)<=0 && x(i)>=-1
                    y(i) = d1f2(x(i));
                elseif x(i)>=0 && x(i)<=1
                    y(i) = d1f2(-x(i));
                elseif x(i)<=2 && x(i)>=1
                    y(i) = d1f1(-x(i));
                end
            end
        case 2
            for i=1:length(x)
                if abs(x(i))>2
                    y(i) = 0;
                elseif x(i)<=-1 && x(i)>=-2
                    y(i) = d2f1(x(i));
                elseif x(i)<=0 && x(i)>=-1
                    y(i) = d2f2(x(i));
                elseif x(i)>=0 && x(i)<=1
                    y(i) = d2f2(-x(i));
                elseif x(i)<=2 && x(i)>=1
                    y(i) = d2f1(-x(i));
                end
            end
    end
end