%
% Example 3 in page 149
% Use Algorithm 3.4
%
n       = 20;
x       = zeros(n+1,1);
fun_val = zeros(n+1,1);
h       = zeros(n,  1);
A       = eye  (n+1,n+1);
rhs     = zeros(n+1,1);

x(1)    = 0.9;
x(2)    = 1.3;
x(3)    = 1.9;
x(4)    = 2.1;
x(5)    = 2.6;
x(6)    = 3.0;
x(7)    = 3.9;
x(8)    = 4.4;
x(9)    = 4.7;
x(10)   = 5.0;
x(11)   = 6.0;
x(12)   = 7.0;
x(13)   = 8.0;
x(14)   = 9.2;
x(15)   = 10.5;
x(16)   = 11.3;
x(17)   = 11.6;
x(18)   = 12.0;
x(19)   = 12.6;
x(20)   = 13.0;
x(21)   = 13.3;

fun_val(1 ) = 1.3;
fun_val(2 ) = 1.5;
fun_val(3 ) = 1.85;
fun_val(4 ) = 2.1;
fun_val(5 ) = 2.6;
fun_val(6 ) = 2.7;
fun_val(7 ) = 2.4;
fun_val(8 ) = 2.15;
fun_val(9 ) = 2.05;
fun_val(10) = 2.1;
fun_val(11) = 2.25;
fun_val(12) = 2.3;
fun_val(13) = 2.25;
fun_val(14) = 1.95;
fun_val(15) = 1.4;
fun_val(16) = 0.9;
fun_val(17) = 0.7;
fun_val(18) = 0.6;
fun_val(19) = 0.5;
fun_val(20) = 0.4;
fun_val(21) = 0.25;

for ii = 1:n
    h(ii) = x(ii+1) - x(ii);
end

for ii  = 1:n-1
    A(ii+1,ii+1) = ( h(ii,1) + h(ii+1,1) ) * 2;
    A(ii+1,ii  ) = h(ii,1);
    A(ii+1,ii+2) = h(ii+1,1);
    rhs(ii+1,1)  = 3 * ( fun_val(ii+2,1) - fun_val(ii+1,1) ) / h(ii+1,1);
    rhs(ii+1,1)  = rhs(ii+1,1) - 3 * ( fun_val(ii+1,1) - fun_val(ii,1) ) / h(ii,1);
end

c = A \ rhs;

b = zeros(n,1);
d = zeros(n,1);
for ii = 1:n
    b(ii,1) = ( fun_val(ii+1,1) - fun_val(ii,1) ) / h(ii,1);
    b(ii,1) = b(ii,1) - h(ii,1) * ( c(ii+1,1) + 2 * c(ii,1) ) / 3;
    d(ii,1) = ( c(ii+1,1) - c(ii,1) ) / ( 3 * h(ii,1) );
end

%
% Plot figure of 
%      S(x) = fun_val(j)+b(j)(x-x(j))+c(j)(x-x(j))^2+d(j)(x-x(j))^3
% on [ x(j), x(j+1) ] for j = 1, 2, ..., n
%
m       = 10;
sub_x   = zeros(n*m,1);
y       = zeros(n*m,1);
Lag_val = zeros(n*m,1);
for ii = 1:n
    sub_x((ii-1)*m+1:ii*m,1) = linspace(x(ii,1), x(ii+1,1), m);
    y((ii-1)*m+1:ii*m,1) = fun_val(ii) * ones(m,1) + ...
                 b(ii) * ( sub_x((ii-1)*m+1:ii*m,1) - x(ii) * ones(m,1) );
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        c(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^2;
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        d(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^3;
    %
    for jj = (ii-1)*m+1:ii*m
        Lag_val(jj,1) = neville(sub_x(jj,1), x, fun_val);
    end
end
figure(1)
plot(sub_x,y,'LineWidth',2)
hold on
plot(sub_x,Lag_val,'r-','LineWidth',2)
legend('Cubic','Neville','Location','SouthEast')
axis([0 14 -10 4])
hold off

figure(2)
plot(sub_x,Lag_val,'LineWidth',2)
axis([0 14 -10 6])