%% Question 4.8.6 page 212

n = 4 ; %(nodes-1)
v = 4*ones(1,n-1) ;
A = diag(v) + diag(ones(1,n-2),-1) + diag(ones(1,n-2),1) ;

f = @(x)(1./x) ;
xn = [1/2 , 5/8 , 3/4 , 7/8 , 1] ;
b = [f(xn(2))-(1/6)*f(xn(1)) ; f(xn(3:n-1)) ; f(xn(n))-(1/6)*f(xn(n+1))] ;
c = A\b ;
c = [(1/6)*f(xn(1)) ; c ; (1/6)*f(xn(n+1))] ;
c = [2*c(1)-c(2) ; c ; 2*c(n+1)-c(n)] ;

h = 1/8 ;
xn = [3/8 , xn , 9/8] ;
x = [1/2 : 0.01 : 1] ;
cTimesBi = zeros(n+3,length(x)) ;

px = zeros(1, length(x)) ;

for i = 1 : n+3
cTimesBi(i,:)= c(i)*exemplarCubicBSpline((x - xn(i))/h) ;
% px = px + c(i)*Bi(i,:)
end

px = sum(cTimesBi , 1) ;

figure(1) ;
plot(x , px , 'r' , x , f(x) , 'k') ;

%% Question 4.8.10 Page 213
n = 5 ; %(nodes-1)
v = 4*ones(1,n-1) ;
A = diag(v) + diag(ones(1,n-2),-1) + diag(ones(1,n-2),1) ; % Eqn 4.49

f = [1.0000000000 ; 0.9181687424 ; 0.8872638175 ; 0.8935153493 ;...
    0.9313837710 ; 1.0000000000] ;
xn = [1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0] ;
b = [f(2)-(1/6)*f(1) ; f(3:n-1) ; f(n)-(1/6)*f(n+1)] ;
c = A\b ;% Eqn 4.49
c = [(1/6)*f(xn(1)) ; c ; (1/6)*f(xn(n+1))] ;
c = [2*c(1)-c(2) ; c ; 2*c(n+1)-c(n)] ;

h = 0.2 ;
xn = [0.8 , xn , 2.2] ;
x = [1.0 : 0.01 : 2.0] ;
cTimesBi = zeros(n+3,length(x)) ;

px = zeros(1, length(x)) ;

for i = 1 : n+3
cTimesBi(i,:)= c(i)*exemplarCubicBSpline((x - xn(i))/h) ;
% px = px + c(i)*Bi(i,:)
end

px = sum(cTimesBi , 1) ;
figure(2) ;
plot(x , px , 'r') ;
hold on ;

x = [1.1 , 1.3 , 1.5 , 1.7 , 1.9] ;
cTimesBi = zeros(n+3,length(x)) ;

px = zeros(1, length(x)) ;
for i = 1 : n+3
cTimesBi(i,:)= c(i)*exemplarCubicBSpline((x - xn(i))./h) ;
% px = px + c(i)*Bi(i,:)
end


px = sum(cTimesBi , 1) 
plot(x , px , 'g+') ;
hold on ;

%% comparison between the values calculated in line 71 above and the actual
%% values from the table 4.8
fx = [0.9513507699 0.8974706963 0.8862269255 0.9086387329 0.9617658319] ;
plot(x , fx , 'k*') ;
Norm2Difference = norm(abs(fx-px))